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SUMMARY 


The present paper explains and illustrates the use of the 
familiar Stodola method of calculating natural modes of vibration 
as applied to coupled bending-torsion modes. At the present 
time, two methods of computing coupled modes are in common 
use, these being the Rayleigh-Ritz method in which the un- 
coupled modes are first computed and then used as coordinates, 
with which the coupled modes may be approximated, and the 
influence coefficient method. The advantage of the Stodola 
method over the Rayleigh-Ritz method lies in the fact that it 
involves only slightly more labor than would be required to 
compute the uncoupled modes by the same method and also in 
the fact that the exact mode shape is found rather than an ap- 
proximation of undetermined accuracy. Further, the necessity 
of solving the determinantal equation is done away with so that 
any additional work over that required to compute the un- 
coupled modes is more than made up for. The advantage over 
the influence coefficient method is that the work can be per- 
formed by a technically unskilled operator, whereas the influence 
coefficient method requires knowledge of matrix algebra. Also, 
the resulting matrix will be of order 2” for m stations on the 
wing, so that for large wings it may be necessary to employ 
mechanical means to evaluate the results. 


SYMBOLS 


vertical displacement of elastic axis, ft. 
mass per unit area, lb.sec.?ft.~3 
spanwise coordinate, ft. 
chordwise coordinate, ft. 
vertical coordinate, ft. 
Young’s modulus of elasticity, lb.ft.~? 
vertical shear force, Ib. 
shear modulus, Ib.ft.~? 
area moment of inertia about neutral axis, ft.‘ 
weight moment of inertia per unit span about elastic 
axis, lb.ft. 
polar moment of inertia about elastic axis, ft.‘ 
= semispan, ft. 
bending moment, Ib.ft. 
moment about elastic axis, Ib.ft. 
moment about centerline, lb.ft 
= vertical load per unit span, Ib.ft.~? 
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torque per unit span, Ib. 

static unbalance per unit span, Ib. 

torque, lb.ft. 

weight per unit span, lb.ft.~! 

angular displacement about elastic axis, rad. 


frequency of vibration, sec.~! 


(1) OUTLINE OF THE GENERAL PROCEDURE 


v | ViiE APPROPRIATE DIFFERENTIAL EQUATIONS (Ap- 
pendix 1) are 


22 Ph] «? 
E& 4 ae 
g 


dx? dx? 


d g 4 m 
dx” dx} 


where x is a spanwise coordinate, / is the displacement 
of the elastic axis of the wing in bending, 6 is the dis- 
placement in torsion about the elastic axis, wu is the 
weight of the wing per unit span, S is the static moment 
per unit span about the elastic axis, J» is the weight 
moment of inertia about the elastic axis, EJ is the 
bending stiffness, G/J is the torsional stiffness, and w 
is the frequency of vibration. 

Eqs. (1) may be integrated in the Stodola manner 
according to the following scheme: 

(a) Substitute into the right members of these equa- 
tions assumed deflection for 4 and J. 

(b) Integrate Eq. (la) twice and Eq. (1b) once in 
the manner prescribed by the boundary conditions of 
the problem (see next section), thus reducing these 
equations to the form: 


(la) 


(1b) 


—"' (Sh + 19) 
g 


EI(d*h/dx*) = M 
GJ(d0/dx) i 
(c) Divide the resulting equations by EJ and GJ, 
respectively, and integrate again, thus arriving at new 
expressions for the deflections. 
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(d) Repeat the process using the computed de- 
flections (each reduced by a common factor if desired), 
continuing until there is no change in the shape of the 
assumed and computed deflection curves. When this 
point has been reached, the constant of proportionality 
between the assumed and computed curves determines 
the frequency. 


(2) BounparRy CONDITIONS 


The modes of vibration of an aircraft wing in flight 
may be treated essentially as the modes of vibration of 
a free-free beam. The ground modes may evidently be 
approximated by those of a cantilever beam. Also, 
because of symmetry, it is sufficient to consider only 
half the wing, say the right half, in either case. 

Taking the origin of coordinates on the centerline 
and designating the semispan by L, the boundary con- 
ditions applying at the tip are 
é d*h = EI d*h -GJ dé = 


— fi — (2) 
dx dx? dx? dx 


x=L, 
corresponding to the physical requirements of zero 
shear, moment, and torque ata free end. These condi- 
tions may be satisfied by performing the first two inte- 
grations of Eq. (la) and the first integration of Eq. 
(1b) from the tip inboard, commencing in each case 
with a value of zero when x = L. 
At the centerline, the boundary conditions to be ap- 
plied vary with the mode of vibration required. In 
the case of cantilever motion they are 


x=0, dhk/dx =h=0=0 (3) 


which are easily fulfilled by performing the remaining 
integrations from the center outward, with values of 
the integrals equal to zero at x = 0. In the free-free 
case, it is convenient to distinguish the conditions of 
axial symmetry with respect to the centerline (sym- 
metric modes) and those of central symmetry (anti- 
symmetric modes). These two sub-cases are defined 
by the relations 


oe = h(—x) 


6(x) = 0(—x) (4) 
jh(x) = —h(—x) 
\0(x) = —0(—x) 


respectively. The boundary conditions for the free- 
free modes are obtained from the well-known conditions 
that, for equilibrium of any free body, the sum of all 
forces and moments must vanish. For symmetric 
modes, this may be reduced to requiring that the sum 
of all vertical forces, F;, and the sum of all moments, 
My;, about the elastic axis be zero, since the condition 
of symmetry ensures the vanishing of moments about 
the centerline. Similarly, in the antisymmetric case, 
the only requirement is that the sum of all moments, 
M,,, about the centerline vanish, since symmetry about 
the origin implies that the sum of all forces and of all 
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moments about the elastic axis will be zero. In addi- 
tion, there are the conditions of zero bending slope at 
the center for symmetric modes and zero displacement 
at the center for antisymmetric, making the required 
three boundary conditions for each case: 


Symmetric: 
. ih 

x = 0, = = 0, UF, = 0, Mo: = 0 (5a) 
dx 


Antisymmetric: 


kh=0, 06=0, )Ma=0 (5b) 


Since the motion is simple harmonic, the forces and 
moments are proportional to the displacements, so 
that the total vertical force over the semiwing becomes 


L 
wf (uh + S0)dx 
0 


the total moment about the elastic axis becomes 


L 
wf (Sh + Ie0)dx 
0 


and the total moment about the centerline becomes 


L 
wf (uh + S0)x dx 
0 


The boundary conditions for the two types of free-free 
motion are thus summarized by the equations: 
Symmetric: 


) 3 
ine Sen ri (uh + S0)dx = 0, 
0 
a 


x= 0, 


dx 
f (Sh + Ipp)dx = 0 (6a) 
0 


Antisymmetric: 
t 
z=0, 4=0, @= 0, [ (uh + S0)x dx = 0 (6b) 
0 


The boundary conditions applicable to a wing carry- 
ing a large concentrated load at the tip such as fuel 
tanks or ram-jet units may be obtained by requiring 
that the shear and moment of the tip load shall be 
balanced out by the corresponding forces and moments 
over the remainder of the wing. This may be ac- 
complished by taking the initial values of the shear and 
moment at the tip to be those due to the tip load and 
requiring that the total force and moment over the en- 
tire wing be zero. 


‘ 


Clearly, boundary conditions of the “integral” type 
just discussed cannot be satisfied by performing the in- 
tegrations in any specified manner. Instead, it is nec- 
essary to perform the integrations from some arbi- 
trarily selected point (say x = 0) and determine the 
constants of integration in such a way as to satisfy the 
boundary conditions. Specifically, in the symmetric 
case, the three final integrations result in deflections 
h, and @,, which differ from the true deflection / and @ 
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by two as vet undetermined constants ao and bo: 


h h, + aol 
6 = 6, + dol 


and these constants must be such that 


, z \ 
[ (uh + S0)dx ff (uh, + S0,)\dx + 
0 0 
L L 
ao f be dx oe wf S dx 
0 0 


oL L 
/ (Sh + [0)dx | (Sh, + [90,)\dx + 
0 


mL - 
a f S dx + wf Ie dx = 0 J 
0 0 


Solving the resulting equations simultaneously results 


II 


rv 


in the expressions 


[nae f" (uh, + S0,)dx is— fr six ff (Sh, sail 
f was fr Ip dx -| fo sax | 
(9) 


, (Sh, + [¢0,)dx =~ f sax fo (uh, + S0,)dx 
0 0 
Fs ‘nae Igdx — ve sasf 
0 


for do and do. 

The problem is simpler in the antisymmetric case, 
since only one constant must be found in this manner 
that which determines the slope of the 
Thus, the bending com- 


—namely, 
bending curve at the center. 
ponent may be written 


h = h, + aox (10) 


Applying the condition 


E 
f (uh + S0)x dx = 0 
0 
L 
f (uh, + S0,)x dx 
0 Se ee (1 1) 


itis seen that 


| lie ra 


px? dx 
0 


(3) COMPUTATION OF HIGHER MODES 


Modes above the fundamental may be calculated 
by the present method provided that, following each 
approximation, the computed deflections are made 
orthogonal to the first mode. This, of course, requires 
that the first mode be computed to a reasonable degree 
of accuracy. 

The condition of orthogonality between any two 
modes, say the ith and jth modes, may be expressed 
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by the formula 

L 
f [uh shy + S(hO, + O:h;) + [98,0;\dx = 0, 1#7Z (12) 

0 
Physically, the above relation states that, if each ele- 
ment of mass (per unit area) of the beam be multiplied 
by its total vertical deflection in space in any two of the 
normal modes of the wing and the product summed over 
the entire wing area, the resulting summation will be 
zero. The relation is thus analogous to the familiar 
condition of orthogonality of the uncoupled modes 
[cf. Appendix (2) ]. 

In order that the Stodola process shall converge to the 
second mode rather than to the first, it is therefore 
necessary that the computed deflections satisfy Eq. 
(12), as well as Eqs. (6). This is accomplished by 
computing deflections satisfying Eqs. (6) as was done 
for the first mode and, in addition, the quantities 


L L 
f h(uh; + S0,)dx; f 0(Sh, + [ef)dx (13) 
0 0 


where /: and @ satisfy the boundary conditions of the 
required mode but not the conditions of orthogonality 
with the primary mode. Having computed the quan- 
tities defined by Eq. (13), it is then only necessary to 
subtract from the deflections h and @ the quantities 
ah; and a,6,, where 


x [h(uhy oo S6;) )dx at A( Sh, + I9;) |\dx 
me ————- (14) 


—— 
yi [ah,? + + 2Sh)6, + 16,7 ]dx 
0 


Continuing in this manner, all the higher modes may be 
found. 


aq = 


(4) DETAILED PROCEDURE AND NUMERICAL EXAMPLES 


The examples that follow illustrate the present 
method as applied to the computation of the first two 
symmetric and antisymmetric modes of an aircraft 
wing carrying tip fuel tanks. The wing has been di- 
vided into eight sections of nearly equal width, Ax. 
The running mass of the wing, uw, has been converted 
into a series of effective concentrated masses by multi- 
plying the value of u at each mid-section by the width 
of the section. In addition to these masses, there is a 
concentrated mass at the center due to the fuselage and 
one at the tip due to the fuel tank. All masses with 
the exception of the last two mentioned are considered 
to act at the middle of their respective sections. The 
unbalance, , and the mass moment of inertia, J, are 
similarly tabulated. The stiffness parameters are listed 
in the form Ax/E/ and Ax/GJ, where the values of 
EI and GJ are taken at the midpoint of each section. 
In this way two operations in the integration process 
are combined. For the antisymmetric modes, the 
spanwise distance to the midpoint of each section and 
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ANTISYMMETRIC ON. 
2 2)/81@| @ 2) a , (2) io). 
Ax |(4Ax} (5Ax) (7,Ax)| 2F (Z/t-Ar(572,Ax) 
Llah3)ax~£[(0)x(2)] 
Joble 
Basic Parameters 
(/) (43) 
del 2, 8, ahrI5a)Ax|(%, = 4 Ax 
J (aki + 25h 64+LG)axr=Sfas)xh + 4x6] 
Table la 
Add i/tone/ constant data far second modes; 
Jymmertt/e or catisymmnerric. 
O10) 210) M1) 77| | @ 
He) A | @ \(nxsiax@hs (oy Ml 5.a\ (edz) 3/2) ial 
(/) | @0)\ (2) | @2)| (23) | 24) | (5) 
M2940) E(20)| Car \2(2) 23) (2Aav 
JTable 2 


fntegration of the bending equation 
Arrows indicate direction in which columns aresummed 


(27) | ror (29) | G02) |G) 1 G2) | 
(2S) shor (28)ay |(29x(7)| E30) | (BAe 


Hl | 


Tab/e 3 
INTEGRATION OF THE TORSION EQUATION 


Oye 


— 


(26) | 
Mx(A) 




















— 





to the tip is needed, which is listed as iy = x L,* and, 
in addition, the quantities uwhpAx and Sho Ax appearing 
in Eq. (11). Since all of these parameters remain un- 
changed from one iteration to the next, they may belisted 
on a separate card or sheet for repeated use. 





*See Appendix (3), footnote * 
tion of this notation. 

+ As a precautionary measure, 
basic parameters be expressed in their proper units. 
any constant factor may be left out, provided the same factor is 
omitted in both equations and provided cognizance is taken of 
the fact in determining the frequency. any 
factor from one equation alone will result in incorrect values for 
the relative amplitudes of the computed deflections. 


(page 262), for the justifica- 


it is advisable that all of the 
However, 
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From these basic parameters, the quantities {pu dx 


Soe Sdx, So" Io dx, 


71 yt L L 

| dx wh wax f Ig dx — | f S ax | ri 
0 / 0 

f sax /$ [lade Ip dx — | f Sax] 

/) ) 

fu afi “ ax f b dx — | f sax 

0 \ 0 0 0 f 


(9) and required for symmetric modes 
Similarly, the quantity 


appearing in Eqs. 
only, may be computed. 


L L 
f uho® dx -{ u(x/L)? dx 
0 0 


which is needed in the antisymmetric case may be 
evaluated. The above data are listed as shown in 
Table 1 

Having listed all the constant data for the problem, 
it is then necessary to assume initial shapes for the de- 
flections # and @. While these may be entirely arbi- 
trary to begin with, it is advisable that they be ad- 
justed so as to satisfy Eqs. (6a) or (6b), according as 
the mode to be computed is symmetric or antisymme- 
tric. This may be accomplished by selecting two en- 
tirely arbitrary deflections, say h, and 6,, and substi- 
tuting them in Eqs. (9) or (11). The assumed deflec- 
tions are listed in Columns (11) and (12), respectively, 
of Table 2. Columns (13) 
by multiplying the assumed deflections at each station 
by the appropriate mass and unbalance, and the shear 
force F (Column 15) obtained by summing Columns 
(13) and (14), starting at the tip with the initial shear 
due to the tip weight. This column is then averaged in 
Eq. (16) for use in the next integration, which is built 
up from the average of Column (15) multiplied by the 
appropriate Ax (Column 2) and summed from the tip 
inboard to yield the bending moment J/ (Column 18). 
This column may, in general, be taken as initially zero 
at the tip, since the moment due to even an extremely 


and (14) are then computed 


large tip weight will be small in comparison to the total 
bending moment. Multiplying the average bending 
moment (Column 19) by Ax/EI (Column 6) and per- 
forming two additional integrations gives the computed 
bending deflection curve, apart from certain constants 
of integrations to be found later in the process. The 
last two integrations are performed from the center out- 
ward and thus satisfy at least one of the required 
boundary conditions at the center in either mode. 


The computation of the torsional component is out- 
lined in Table 3, and proceeds along similar lines. 
Columns (26) and (27) are computed from the assumed 
deflections and the given unbalance and inertia param- 
eters, listed in Columns (4) and (5). These are sum- 
med starting with the torque at the tip due to the tip 
weight and proceeding inward. The resulting torques 
are then averaged and multiplied by Ax/GJ [Column 
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— TED DEFLECTIONS 
BENO/NG TORSION 
“| | (33) G34) | (35) G6) (37) (8) 

| ( (25)x(3) | (25)x(@) | o2.(2)| @2x()| Gs) *, |cox)+6, 














X03) 36d S(5)= Xe)= " 
a,* ad 
| 
“ Table 


CALCULATION OF THE FINAL /NTEGRATION CONSTANTS; SYMMETRIC 700ES 


the result being integrated from the center out- 





i) iy 
ward in Column (31). 

The final boundary conditions of the problem may now 
For symmetric modes these are determined 
from Eqs. (6) and require that Columns (33), (34), 
(35), and (36) be added to the schedule listing the quan- 
tities uh, Ax, Sh, Ax, SO, Ax, 190, Ax. These quantities 
are computed from the basic mass parameters multi- 
plied by the deflections h, and 6, at the stations where 
the masses it is therefore 
necessary to know not only the deflections at the center 
and tip but also the deflection at each mid-section. The 
required deflections are listed in Columns (25) and 

(32). Summing Columns (33) through (36) and sub- 
stituting the results in Eqs. (9) gives the required 
values for dp and bo, and these in turn are substituted 
in Eq. (7) to determine the final deflections satisfying 
all the boundary conditions of the symmetric mode 
(Table 4). 

The final boundary conditions for the antisymmetric 
mode are simpler and require only the addition of 
Columns (39) and (40) listing the quantities uh, Ax 
and s6,h9 4x, found in exactly the same manner as the 
corresponding quantities of the symmetric mode. The 
summation of Columns (39) and (40) are then sub- 
stituted into Eq. (11) to obtain ao, and this in turn is 
substituted in Eq. (10) to yield the final bending deflec- 
tion. No constant need be determined for the tor- 
sional component of the deflection (Table 5). 

When the first mode has been computed to a reason- 
able degree of accuracy, the second mode may be found 
by augmenting Table with the data shown in Table 
la which were easily obtained from the first mode. 
(13) and 
flections may be obtained satisfying the boundary con- 
ditions (6), as well as the condition of orthogonality 
with the first mode. These deflections are entered in 
Columns (11) and (12), and the process is repeated 
exactly as though the first mode were being computed, 
except that the computed deflections must again be 
substituted into Eqs. (13) and (14) to obtain the quan- 
tity a, and the final deflections. This requires the ad- 
dition of the columns in Tables 4a and 5a to Tables 4 
and 5, respectively. 


be applied. 


are assumed to be located; 


Then, by means of Eqs. (14), assumed de- 


After several approximations, the coefficient a, be- 
comes increasingly negligible. 
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Table 4a 
REMOVAL OF THE FUNDAMENTAL MODE-SYMMETRIC CASE 























COMPUTED DEFLECTIONS 
BLENDING JORS/ON 
i (39) (40) | @/) 42) (32) | 
Lal os. (9) | G2)x(/0)| 2.6) | @5+E/) 
539- X40 
a,= 
Table 5 


CALCULATION OF THE FINAL INTEGRATION CONSTANTS; 
ANTISYMMETRIC Moors 
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t(45)-  3(e)- 
a, = 
Table 52a 
REMOVAL OF THE FPUNDAMLNTAL MODE -ANTISYMMETRIC CASE. 
(SYMMETR/C) (ANTISYMMETR 7 
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oy Ax fdednid Axx oO (gx) <ib Be A > 
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(8) |/.66667| 0/4 | .0/03 | .020 \/50000|2.38/0 | 25763) .0/34/ | .00996 
Sea & oi | suena 
7/P|\—— |4./86 | .559 | 7757. — /.00000Y. 18600 | .55900 
POTALSV 2.66667) 7.3// paace7ser2.736) £ ahsdx = /.58682x 10" 











dx = 3. as eo - 3. > 2/2? io? 
[adx = 7.3/1 x10"; [Sex =-306275 x10"; [Jax =2/2.736%/10 


(E “adx) x/0* 
| Faax/', dx-[[‘sax] 
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Frucx [ax] saxq 


29 








= 0047: 


-00/93/ 


‘wg dx) x/o” 


=./376/ 
fadrf lax [fsax) 7 
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| * units for 44x are pounds, SAx pound feet, Lx pound feet squared. 
| Ax ss infeet, Eland G/ in pound feet Squared; the foctor J" /acking 

| fa the mass parameters /s supp//ed when rhe frequency /§ found. 
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| (a) 70/023 | 5/534. 25/239 4, °3|40.6/ 038 | 85.43608 /G5./3009 
ww feet Eee Cee Ce De od 
First Symmetric A 
Wa [ Ge) Bod [ eo) | 2) | Go | Ga (32)¢ 
ye () )x(4) | | V2)x(5) [S[iedzd]| (20)av | (29)x(7)| = (30) | (s/av 
aie 2688%ro0e0e —— || —|_@ | o_| 
Lf) |-. 07632 |- 0/368) 664/2 se ‘| 02394|.023 4 | 0/197 
(2) |-05/66 |- /770| 1: | -78880| 09654) ./2048| 07221 
Fa fase [0 2/5701 8 02340] .06939) 09/6 | 30962| .2/506 
} |-.0/264 |-.00744| 87 528 | 88532 | 36892 | .66956| 48970 
(5) | .o0ase 2 00/98 | .. 09536 | .09409| 738 17,325.90 99723 
6) eens eee 89282 | .87730 |/ 0|2.51130 |/.9/ 860 
(7) | .02376| .00278| .86/78| .9435/ |2.2/970 |4.94/00| 3626/5) 
@) | .00886! onan 32522). B205/ |7.95364|6,694645, 7/762| 
| 7/2 | 25900 | 125678) | 8/578 | ——| i ms le.cose¢ 
Jab/e 7. 
Integration of the Torsion Lauation 





« The, assumed deflectk are free of the s¢r0 mode; 
Aence the 1 athad adjust mene. mentioned inthe 3-d par agraph, 
Section 4,75 pot NeCCEsSSAY. 


“x The values in these spaces should be opproximately egual 
fo the negative of fhe shearand torgue af the hase station. 
da this meamer the correctness of the values of the constants 
a and b/s verified before proceeding fu riper 


+ Values in oar oe and (32 @re avera re Asa VecHons, 
with the exceptions of hose atthe 5asean Le tip wich 
are the actual deflections at these s/auvo aS, 


When the computed deflections and the assumed 
deflections differ only by a constant factor, the process 
has converged to the desired mode. By Eggs. (1), this 
factor determines the frequency of the mode, provided 
account is taken of constant factors omitted from the 
basic parameters. A convenient check as to the con- 
vergence consists of reducing the computed amplitudes 
by a constant factor such that at one station the ampli- 
tude of either the bending or torsion component is 
identical to the corresponding assumed amplitude. 
At convergence then, all of the reduced amplitudes of 
the computed deflections will be equal to the corre- 
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First Symimet. Yode 
es = = 2 = === os of = == 

| \conmsure 0 DEFLECTIONS | NORMALIZED" DEFLECTION 
| | aenoiwe | TORSION BENDING TORSION 

| (33) (34) 35) | (36) (37) (38) A, a, 

cr. y ao , | o / f ) sf f- 5 ‘4 

Wa \(25)x(3) | (25)x(4) | (32x(4) | G2x(S) | (25) + |(B2) +, | En LENOLE 
Asse i) 7) 7) oO 38/2232 |-/.42/10 \ - 25873} L0%4676) 
—tr—- ++ : — 

‘/) pad 981 | 23/75 20%) 07) /2 P. i= 35237 (409/35 ~ 25 356| ~. 009566 
(2) | 2 353 74, .0524! | 0/603 (3958 | 34 20262 '- 34889 | ~.-23272 "009/57 
head ed Es tases ae 
1 (3) £03907 | 499173 .Ab/@3| 4, v2e7 | -27. 7.73249 | 20604) - = 008/87 | 


| (4) | 9.22743| 2.6364/  .05770| .576; a | 5.77987) ~.93200| ~./07/2 | -.006327 
1 (5) | 8.4/075| 678778 ./60/6 .70¢08 4/4278 | - .42387| .028/2\-.002877 
> + —-— + t 


(6) 10.65720' 3.03289 23580 8446/8 13537560 | 43750\ .280/4\ 003377 
(7) | 4 27386 | 7.32775 | .2239/ 67084 |. | ea.sesv2 2.20505 | .54676| .0/4%8 





(Sehee eh A ne cabeteeisiaah | | 
\(8) | 2.3//83| 470005 [05009] "836 (2100857 | 4.29672| 8662/6) .029%7 


oars - ™ T — —4 
71 b/9.927/9 \103.65887 | 3. 74230 |97.07337/ 473/376 | 5.27354 00000) 035738 
FOTASED, 2. 2048)/ 84.22 2038 | 4.53962 p62 |57. 35506] 

















("th, dx = 269.92088; Z “SQ, dx = 4.53962; =/4.78 Sec 
“o 





rhe paw 
4 >> a: ’ = 4 586 
J, 44, Wx = /34.22034 ; L, 9. AX = 51,3558 


&, = ~ (./376 1269.8 20438 + £,53962) + (-.00/98/)/34.22038 + $1, 35SEC) = »38,/2232 





5, =- (.004729)(/34. 22034 + 5/.35$86) + (-.00/98/)(269.82048 +4.5 3962) #/42//0 2 








7ad/e 7c 



























































1, a3) | a) 
Ne h, 8 (22,75 AJAX Sh, +4, G)Ax 

PASE | ~.25878 |-.00964676) ~ 82207 | 66 3// | 
/ ||—.25356|-.002566|-.20700 |-.09000 | 
2 —23272 ~ 009/57 —./4L692 e 06936 | 
3 |-78825 |008787 | - .09287 F-. S179 | 
4 |-./07/2 |-.006327|-.04499 |-.020/2 
S || .028/2 |.002877| .005/4 | .0 10248 
6 || .240/4| .003377| .03524 .03099 
7 | .see7e| .oa96e| .0szee | .03653 
8 || .862/6| .029/67| .0/237 | .00946 

TIP | 400000) 035798 | 4.2060/ | .8/499 

ys “(ah7+25h, Gt1, 6; ) dx =/.6/3 25 








Table 8 
Additional constant data required 
for computation of Second syamerric more 


sponding assumed amplitudes. This process is some- 
times referred to as ‘‘normalizing.”’ 

In the examples shown in Tables 7a-7c, 8, %a, 9b, 
10, 11, 12a, and 12b, the deflections have been nor- 
malized so that the bending deflection is unity at the 
tip station, with the exception of the first antisymmetric 
mode where it was elected to make the tip amplitude 
of the torsion component unity because of the relatively 
small amplitudes of the bending component. 

The computation of the first symmetric mode 
carried out in detail showing all of the necessary steps, 
while the remaining modes are exemplified by showing 
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iy (//) (/2) (25) | ©2) | (33) | G4 | ©) | Ge) 
we | 2 8 C4), | GAav_| (25)x(3) | (29x) | (52)x@)| (32)x) 
BASE | ~.04203 | 05153 O O O O 7) o 
(/) |-.0394 | .0483 24/3 | ~./658 /9245\ .07263 | -,0499/| ~-.23709 
(2) || -.02// . 0329 1/273 | -.9726 735394| .26580| —.2/592 | -/.58004 
(3) | .0/7% | -.00/7 32/26 | -2.7832 | /558//\ €/§86 | ~53354\ -$.33539 
(4) } .0865 | -0658 6.825/ |-6,/34/ | 2.8/877| .60536| ~.72382 | -72505/ 
(5) | ./998 | -./804 /2.7A7a \J2./340 | 2.53673 | 2.04723 | ~/.94872 | -8.59087 
(6) | .3753 | -386/ 2/.2233 |-22.8966| 3./7888 | 2.69437 | ~2,81399 | -/0,07A50 
(7) | -6277 | -.7745 35/183 \-A3.214/ | 3,33624| 2/6856 | -2.66847| -7.9946/ 
(8) | .6868 |-4260/ 48,665/ |-68.6/64| .68/3/| .$0/25| -.70675\ -/.37233 
7/P || 40000 |-/.4200 545824 |\-SQCL52 |€4,73473 |30.51/56 |-A5.08067|-576.69383 
72,77296 |32.67992 |-50,74179 +6/2.429/7 
fh, dx=79.77296 56, dx =-54, 74179 
J Sh, WX = 329.67992 [99.dx=-612482917 
8, = -(./376/)(72.77296 -54 74/79) + (00198!) (39,6 7992 ~-6/9.42911)=-2.29594- 
4, = -(.004722)(39,67992~ 6/9,4Z917) +(-.00/98/)(79.77226 - 54. 14/79) = 2.622/6 























Table Qa . 
Final Boundary Conditions 


only those steps in the process which are unique to the 
particular mode under consideration. All steps omitted 
are carried out in a manner identical to the correspond- 
ing ones involved in computing the first symmetric 
mode. 

The computed mode shapes are shown in Figs. 1 
through 4. 


APPENDIX (1).—DERIVATION OF THE EQUATIONS OF 
MOoTION 


In deriving the differential equation of the beam, 
it is convenient to define a rectangular coordinate 
system so oriented that the x-axis coincides with the 
elastic axis (positive to right), the y-axis coincides 
with a chordline of the beam (positive aft), and the 
axis extends positively downward in a direction per- 
pendicular to the plane of the beam. In addition, 
it is convenient from elastic considerations to define 
two additional coordinates h and @, representing, re- 
spectively, a vertical deflection of the elastic axis and 
an angular displacement about the elastic axis. The 
telation between the two systems is at once seen to be 


z2=h-+ y6 


at any value of x, where / is positive down and @ is 
positive if the rotation takes the positive y-axis into 
the positive z-axis (Fig. Al). 

From the theory of structures,' the deflections / and 
6 due, respectively, to a distributed load per unit span 
P and a distributed torque per unit span Q are given by 
the equations 


(d?/dx*) EI(d*h/dx*) 
(d/dx)GJ(d0/dx) 


| , 
\ (A1:01) 

— 

In the case of a beam vibrating harmonically with a 

frequency w, the force and torque of an element of mass 

per unit area, m, are, respectively, 


OP/Oy = w*mz 
0Q0/dy = w*mzy 


(A1:02) 


and the total force and torque on a spanwise element are 


P= ot fms dy 
Q= ot f mey dy 


in which the integral extends over the spanwise element 


(A1:03) 
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1)| (7) | GA | 5) | Ge) | 7) | Ge) | @) | GO|] Zz E 
574 (25) xa, | (32) 4, | (37)x03)| 6x08) 2k, | 24 |6)-@i|Gd-43| Bes | Le 
BASE| -2.29594| 2.692/6| /.88742| -/.78520| .00762| .00028\-2.30356| 269/88 -.00403 | 051454 
(/) | -205¢64| 25263¢| .4253/| -.22737| 00746| .00028\-2,062/0| 252608) -.03942| .04829 
| (2) | -109864 1.7/95¢| ./5a9¢| -./927| .00685| .00027| -/./0549| /.7/929\ -.021/3 | .03286 
(3) 9c6c| -.09104| -.085/3| .0047/| .00558| .o002z%| .9///2\ -.09/28) .0/742\-.00/74 
(4)| 2.529/6\-3.4494| -.20377| .06925| .003/5| .000/9| 4.526¢0/|~3.442/3| 08657 | -.06580 | 
(5) | 10.45146|-943/84| .06372| —.02332| -.00083| .00008) /0.45229| -943/92\| ./9972 | -./8029 
(G) | /9.6273C}20.20044\ .69/67| ~-.626/4| -.00707| -.000/0 | /9.63443\-20.20434\| .3753/ | —38620 
(7) |32.62236}40.52/94| /.73499| ~/48027| ~.0/609| ~.00044|32.83845+20,52/50)| .62770 | -.7 7456 
(8) |46.369/6+ 6592424 57359| -.62364) -.02537| -00086|46.39453|6592338) .8e682| 7.260// 
7/P | $2.28C46\-77, 95304) 63,05799\-63,53095| -.02943| -.00/05|$2.3/582|-77,95/99| /00000 | -/.49003 
3a.2947c\6a.s0zze| 
Oy = Forex 10* 
[lakh + SRO) dx =68.29476 = 248sec? 
Lo S0h, + J, 06) 2X = -68,34226 
me me ee 
7a@ble_ Ib. 


LReECMOVAL S first Mode 


under consideration. Substituting z = h + y6, and 
noting that # and @ are functions x only, Eqs. (A1:03) 


become 
wn fm dy + af my dy | 
wh f myay + af myray | 


in which {° m dy is the mass per unit span, {© my dy 
is the static moment per unit span, {” my* dy is the 
mass moment of inertia per unit span. The differential 
equations of the motion are, therefore, 


(d?/dx®) EI (d*h/dx?) = w%(uh + S80) 
(d/dx)GJ(d0/dx) = —w®(Sh + Iy6) 


where p = { mdy,S = f mydy, andl, = f{ my’ dy. 


(A1:04) 


t (A1:05) 


APPENDIX (2).—-PROOF OF THE ORTHOGONALITY OF THE 
CouPLED MopDEsS 


To demonstrate the orthogonality of the coupled 
modes requires that it be shown that 


Jf mse y)2;(x, y)dx dy = 


where z, and z, are the total vertical displacements in 


0, i# 7 (A2:01) 


any two modes arising from the same set of boundary 
conditions and where the double integral extends over 
the entire area of the beam. Writing 2; = (4; + y6,), 
zy; = (Ah; + y6;) and integrating with respect to y trans- 
forms this relation into that given by Eq. (12): 


fl 
J [uh sh; + S(h 0; + 6,h;) + I96;\dx = 0, t ca ¥ 
0 
(A2:02) 
It is thus sufficient to demonstrate the validity of Eq. 


(A2:02). If w; and w, are the frequencies associated 
with 2, and z;, then from Eqs. (1), 


= : - = wt(uh, + 50) | 
= Gs “ = —w,°(Sh; + 196,) 

- 
¥. ‘i o acai ills (A2:03) 
. GJ “ = —a,?(Sh; + 118;) 





4 
Furthermore, the boundary conditions associated with 


the various modes discussed may be written in the 
generalized form: 
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“| MM | ~@ (25) | GO| 32 | 4 | & (42) h, 6, 
2 (32 
we | A e (24)4,| Baw | (25)x(9)| 62x00) (8) x2. | 23) G)) BEtzo|_B2Iz- 
BASE| O Oo | 1 oe | @ 7) O oO oO O O 
~ a 
(/) | -.042/| .00¢/4 | ,0097| .2302| .00053| .0087/|-2,3458| -233¢/| -.04/9| .004/4 
} ———~4 
(2) || -/256| O/é2E | 0626} .9070| 00859| 04082 -6.9637 |-6.894/|| -/236| 0/626 
(3) || ~.2026| .03872 2748)\ 2/600\| .04433| /3772\-/, §086|-1/.2338| -.20/4 | 03872 
(4) | ~.27§5| .O80// | .0470| 4.4693| /6305| 2456/\-/6./265 |\752795\\ ~.2732\ 080// 
| 
(5) || -.3342| ./§4/3 2/230 | 8.5985| .25240\| .82503-20.67/4 |-/8,5484\\ -.3325 | ./54/3 
(G) | -.3694, .2870 4 $406\/6.0022| 7708) | M25EB-25. 0676 \-20.5290 || —.3680| .2870 
(7) || -.3779| .5379 GA29C|300070| 68208) /.5787}+29.467A |-2/.0378 || ~.377/ | .£379 
(8) | —.3€72| .A5/5 /2.6285\47,.5037\ J6935| ACOIN=33./328 \-20.5043\| ~.3675\ 85/5 
7/P || ~.3595| 7.0000 /4.4973 55,7870) 7./9380\3/./8493\-34. 5988 |-20/0/5 \\ -.3603| /.0000 
18.99 113 \35.9/086 J 
2 aval 32.2 
[ah h, ax =/8.991/3 = 24,03 sec.! 
oO 
Z 
[ 58,4, IX = 35. 9/086 
Se al 
18.99/13 + 35. 9/086 
ale 7.58682 ~ 54, 5287S 
| Table /0. 
FINA! BOURAAPL CONAAIONS 
| (/) a (43) (ga) Multiplying Eqs. (A2:03), respectively, by /; dx, 
ras &, | @ (dh +5Q)AX(Sh, +1,8)Ax 6, dx, and h, dx, 0, dx and integrating between the limits 
laase| ogre gage ¢ 0 and L, 
i hoe a 
(7) | » | otra) -.0320 | -.0067 Ld dh, i 
7) — h; In? EI ng dx = wef (uh hy + SO:h,j)dx 
| 4 i 3¢ /6LE Cléc | 0040 0 ax 
| y | ca | -aeale Ia) ~ 7L 2 2 L 
LF tere | 05072 | ‘1 y| 5 = a em fi (uh,hy + Sh,b,)dx 
(a) | -.2739 | .aaor | ~.703 0624 0 dah ae 0 
—+ — L 8 
(5) 3325 | /$4/3 | OLI5S | .O557 | f 6, d gy mad dx = - ot f (Sh 0; + [90,0,)dx 
(S) 3 30 | .28 70 | -.0/38/ er | 0 dx dx 0 
create Torso cvcallll It “L L 
(7) | -.377/ | 5379 | -.0026 | .0762 | 0, q iJ ” dx = -o7 ff (S0jh; + 196,0;)dx 
— —>— mem : 0 dx dx 0 
(3) 3¢ / 4 | BS/S 0037 0/32 \2:05) 
7/P -.3603 | 1.0000 | ./3/7 | 6.9496 
| Fs (2h? + 254.8 eax =7 O625C and an obvious integration by parts of the left members 
ye . of Eqs. (A2:05), together with the generalized boundary 
P oo cd iia conditions (A2 :04), yields the relations 
Gadiuofa! COM STA, 4 we QUME 
forcom VR OUIAL ON Of. Se “COW y. dé A eb 


xX 


d*h dh d*h 
(m( 5 BI dx? za) (=) Ea) = 


fx = 0 
(67 a) ae lx = L 


(A2:04) 








Ld dh; > - th 
fue ul 2 & dx - fund ae ax| 
d dé, d dé | 

7J - —— — tJ at x 

f % dx ’ dx - - [ a dx ” 


(A2:06) 
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1h dx, 
limits 
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f_« > 
(@;" — w;*) 


"L 
(w,? — wi) f ph hy; dx = 
0 


ot f Sh, dx — ot f SOsh; dx | 
0 0 


WL 
{w, — wi) f 1900; dx = 
0 


r (A2: 


PL “7 
of SO;h; dx — ot Sh 0; dx 
0 0 : 


°L 


(uh hy -L I900;)dx = 


0 


whence, combining the equations, 


STODOLA PROCESS TO THE COUPLED MODES 


“) | (4) (2) (25) | (32) | 9) | @) | @) | &) 
Jel ok o (24)av | Cav | @5)x(2)\(32)x(0)| 6) x& | (25"E/) 
BASE 0 O 0 0 0 O O 0 
(/) | --298/ | -.002/9 O752| —.0465 004/4.\ —.00025|-6.4062 | -€.3303 
(2) | -—.870/ | -.008/¢6 388) —./735 06C48| -00775\~/9.0/73 |V8B.ATOS 
(3) |-/.3804 | -.0/748 2.//40\ -.3722 »34/03| -.02373 |\-3/.422/ |\-29.3/5/ 
(4) | -4.7666 | -.03075 6.5206| ~6550| /.25522| —03603|-46.0402 \-37.5/96 
(5) | -4.8834 | —.04704 /6.4450\ -/.0034 | /.9§5/1§ | -.09628\-56,4520 £00070 
(SG) | SA~S552/ | -.0626/ 35.5256| /.3377 | 3.73232\ -./9/2 684632 392.9376 
(7) | -.<500 | -.07489 66.66/6|-/.602/ | $.39359| -08425'|-80.4734. \-/3,8//8 
(8) 4753 | -.08357 100.57.96\-/.7882\ /.34877\ -.0/764 \-904833 | 0.0963 
7/P | /.0000 | ~.08696 11S. 7362\ -/.86/3 \/37.263/3 | -/.04075 944867 | 21.2495 
V5 1.35983\ -/ 42650 
F ah fe AX = /5/, 35983 
(-] 
f SE h, ax = -/.42650 
@,=- LSS IOS Rese 50 =~ 94, 48667 
Table /2a. 
FINAL BOLRAAPY CORAlmo/ns 
Thus APPENDIX (3).—PROOF OF THE CONVERGENCE OF THE 


The convergence to the coupled modes of the Stodola 
iteration procedure is dependent on the fact that any 
arbitrary displacement in space of the wing may be 


07) 


the wing. 


follows: 
. . | 

Let the arbitrary displacement be 2(x, y) and assume 

that z can be written as 


2(x, y) = 20(x, vy) + aizi(x, v) + d2%2(x, y) + 


. Oeil ®, ¥) 


expressed as a series of the normal coupled modes of 
This may be demonstrated formally 


as 


(A3:01) 








“L L 
(w;? al ot) f (Sh 0; dx + (w,? —_ wo) f SO;h; dx 
0 0 


(A2:08) 


and since it has been assumed that w; ~ w,, it follows 
that 


L 
; [uh hy os S(hp; ob 6,h;) _ I96 0; \dx => 0, 1 gt j 
) 


Here, 20(x, y) represents the rigid body mode of the 
wing and 2z,(x, y) represents the mth normal mode. 
Multiplying the equation by mz,(x, y) and integrating 
over the area of the wing gives 


ff msc y)in(x, y)dx dy = 
“ aw ff main, y)dx dy (A3:02) 
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(/) (45) (46) (£7) (48) (49) (50) ~ | - | 
SA \(42)x@3) |(s0<(48)| 24, | 34 (e2)-47) (32) GO| a= | GPs 
BASE 0 0 0 Oo 0 0 O o | 
(/) | .20257| .0003/\ 0006 | -.000/ |-6.3309 | —.0464 | -.2980 | -.002/8 | 
(2) | 4334/5| -.00069| 00/8 | -.0007 |-/84803\| -./733 | -.8699 | -.008/6 | 
(3) | 2.647/5| -.0/325| .0030 | —0006 |-29,3/8/| ~.37/6 |-/.380/ | —.0/749 
(4) | 3.80703| -.04087| .0040 | -.00/2 |=37.5236| —.€538|-/.7663 | -—.03078 
(5) | 466029| -.05589| .0042 | -.0023 |-£0,0/12 | -/.00// \-/.8634 | -.047/2 
(SG) | .596/7| -/0849| .0054 | -.0042 |-32.9430| -/.3335 | -/.5507 | -—.06277 
(7) | .0359/ | ~./2208| .0056 | -.0079 |-/3,8/74\| -/,5242 | -—.6504 | —07504 
(8) | .23736| -.0236/| .0054 | -—0/25 | /0.0909|-/.7764| .4750 | -.083G2 
7/P | 2.79856|72.93877| .0053 | ~-.0/47 | 2/.2442| -/.847/|\ /.0000 | — 08695 
13.1939/9 \}-/3.30334. 
; Man 777 kel 
wh (ahh + 548) Bx = /3./99/9 = 39.94 sec? 
Lf (5 Ok, + £09) dx = -/3.30332 
a, =—& A ee = -.0/47 














7@b/eE /ZD 
REMwaA Of SIP?ST Mode 


since all terms other than 2z,(x, y) vanish by virtue of 
the orthogonality. 


Writing the total displacement z in terms of the co- 
ordinates h and @ results in the expression 


h + yO = (aoho + boyOo) + ai(hi + yO) + 
A2(h2z + yO.) + ...dn(An + yO.) (A3:03) 
the coefficients of ho and 6 being different, since no elas- 
tic restraint exists between these motions.* Since 
Eq. (A3:03) is true for all values of y, it follows that 


Il 


Ao + ayhy + aoho +... Anhty\ 
bo + ai91 + a2O2 +... GnOnf 


h 


where do and dp are given by Eqs. (9) or (11), and 





* For cantilever modes, 4) and 6 are both zero. For sym- 
metric free-free modes, 4) and 6) are constants which in the 
present case may be taken as unity. For antisymmetric free- 
free modes, 0 = O and hy = x, or, if ‘‘normalized”’ in the sense of 
the first paragraph on page 264, hy) = x/L, where L is the semi- 
span. 


f ri m2(X, V)Zn(x, y)dx dy 
a, = —.— 
J [mai y)dx dy 


"Lh 
J (uhh, + S(hO, + Oh,) + 1900, |dx 
0 > 
= —— (A3:06) 


j Ss 5 ~ ; 
f (uh? + 2SinO, + I902)dx 
0 


In the case of the Stodola process, the assumed bend- 
ing and torsional deflections may be taken to be of the 
form 


h 


ayhy + adhe +... “ (A3:07) 
0 = a1, + a8, + ... AnO| ~ 
in which the rigid body terms are absent. This is 


because the assumed curves are supposed to satisfy 
the integral boundary conditions of the mode, which is 
equivalent to the deletion of these terms from the series 
representation of the deflections, as may be easily veri- 
fied. Substituting the deflections 4 and @ into the left 
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side of Eqs. (1) and integrating will result in computed 
deflections h“ and 6, which may be written 


, hy he hy, 
so =a,— +4—7 ...&— = 
@)” W,” Wn” 
a . | 
l w,\° w,\" | 
" | av oF ( ) Ashe + ... ( ) aus 
Ww" We Wy 
> (A3:08) 
6 Oo 6, } 
go = a, = > i —— a, — = 
@)" We Wn” | 


l Ww) 2 W@W) 2 | 
= a6; + A0. + ... AnD, 
@1~ We Wy J 


since any one of the individual components of the series, 
say h,, 0,, would yield deflections h,, w, 6,/w2. Con- 
tinuing the process, it is clear that after r approxima- 
tions, the computed deflections will be 


1 an 2r wy} 2r 
h® = — | ah, + Aghg + ... a,hy 
a We Wy 
1 w1 2r w 2r 
ow = zs 110, + 1205 + =" oe Anby 
oF i we Wy, ] 


(A3:05) 


and since w. > a1, @, > w, etc., it follows that the 
terms of the series involving modes higher than the 
first become increasingly small and that eventually the 
computed deflection will approximate the first mode 
to any desired degree of accuracy. 

In a similar manner, it is seen that, if the terms in- 
volving the first mode are removed from the series by 
first determining the coefficient a; and then subtracting 
the terms a,/ and a,6, from the respective components, 
the process will converge to the second mode, etc. 
This ‘‘sweeping’’ must be performed after each itera- 
tion, since inaccuracies or errors in the computation 
tend to magnify the lower modes for the same reason 
that the higher ones are diminished.’ 


REFERENCES 
! Timoshenko, S., Vibration Problems in Engineering, 2nd 
Ed., p. 325, Art. 53, p. 331, Art. 54; D. Van Nostrand Co., Inc., 
New York, 1937. 
2 Beskin, L., and Rosenberg, R. M., Higher Modes of Vibration 


by a Method of Sweeping, Journal of Aeronautical Sciences, Vol. 
13, No. 11, pp. 597-604, November, 1946. 
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U.S. Naval Ordnance Laboratory, White Oak, Silver Spring 19, Md., Attention 

















Elastic Stability of Simply Supported Flat 
Rectangular Plates Under Critical Combina- 
tions of Transverse Compression 
and Longitudinal Bending 


NORMAN GROSSMAN* 


Republic Aviation Corporation 


SUMMARY 


The elastic stability of a simply supported flat rectangular plate 
subjected to combined bending and compression was investigated 
theoretically in an effort to derive the correct interaction curves 
for various plate aspect ratios. By use of the Rayleigh-Ritz 
energy method, stability curves were obtained as a function of the 
plate a/b. The results indicate that the interaction equation 
recommended by ANC-8S is incorrect, yielding optimistic values 
for plates of a/b < 1.05 and conservative values in all other cases 


INTRODUCTION 


ben STABILITY OF PLATES under the action of com- 
bined forces in the plane of the plate for different 
boundary conditions has interested numerous investiga- 
tors. Complete solutions are available for various 
combinations of loads and different boundary condi- 
tions. Among the problems investigated are the in- 
finite strip? and finite plate* under shear and longitu- 
dinal direct stress; infinite strip‘ and finite plate* under 
shear and transverse direct stress; infinite strip’ and 
finite plate® under transverse and longitudinal direct 
stress; finite plate under shear and bending;' finite 
plate under longitudinal bending and longitudinal direct 
stress. ' 

To supplement the above investigations, the elastic 
stability of a rectangular plate under longitudinal 
bending and transverse compression was developed. 
This problem is of importance in the design of box-beam 
webs. 

In the analysis presented herein the Rayleigh-Ritz 
method for computing critical loads was used. Matrix 
iteration was also employed in an effort to evaluate the 
relative effectiveness of the Fourier coefficients em- 
ployed in the Rayleigh-Ritz procedure. 


NOTATION 
a = length of plate in direction of bending stress 
b = width of plate in direction of compressive stress 
2. = wave length in the direction of bending stress 
8 = aspect ratio \/b = a/mb 
t = thickness 


= Poisson's ratio 
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E = Young’s modulus of elasticity 

D_ = plate rigidity in bending, Et*/12(1—v ?) 

x = plate coordinate in direction of bending stress 

y = plate coordinate in direction of compressive stress 

w = plate buckling deflection normal to plane of the plate 

Nz = bending force per unit length acting on plate middle 
plane 

N, = compressive force per unit length acting on plate middle 
plane 

ky = compressive stress coefficient 

ky = critical bending stress coefficient 

y = eigenvalue, r?/16k,8? 

V = internal energy of deformation 

T = external work of applied forces 

dy, = Fourier coefficients 

R_ = stress ratio 

J = unit matrix 

{] = asquare matrix 


Q) = a matrix column 


METHOD OF RAYLEIGH-RITz 


The Rayleigh-Ritz method is usually the simplest and 
most direct procedure for solving elastic stability 
problems with simpiy supported boundary conditions. 
The critical force is determined by minimizing the total 
potential energy in the system. If V is the internal 
strain energy, 7 is the external work and A is the total 


potential energy. 
A=V-T 


so that the first variation 
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Fic. 1. Plate coordinate system. 





the 
rel 
on 
en 
the 
ing 


leis 
wit 
ser 
pot 


Bot 


wh 


ise 


Sin 


Eqs 


(Mn : 


Let 


te 
iddle 


iddle 


ind 
lity 
ns. 
tal 
nal 
tal 





ELASTIC STABILITY OF 
6A = 6V —6T=0 
6V = 6T 


In the analysis it is assumed that the deformations of 
the plate mean plane and the forces on the boundary 
remain invariant at the onset of buckling so that the 
only contribution to 6V comes from bending strain 
energy. Since 6V and 67 can be expressed in terms of 
the deflection w, the problem is reduced to that of find- 
ing a deflection function w that minimizes the functional 
A. 

The unknown function w is constructed by the Ray- 
leigh-Ritz method, expressing w as an infinite series 
with undetermined coefficients. Each term of this 
series must satisfy all the boundary conditions. The 
potential energy is then computed and minimized with 
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respect to the unknown coefficients, resulting in an 
infinite set of linear homogeneous simultaneous equa- 
tions in the coefficients. By setting the determinant 
of this system equal to zero, an infinite polynomial in 
terms of the buckling force is obtained, the smallest 
root being the required solution. 


Since it is impossible to evaluate an infinite deter- 
minant, some of the Fourier coefficients are assumed to 
vanish. Therefore, the critical force obtained in this 
manner represents the solution to a problem where the 
buckled surface is prescribed by a finite set of coeffi- 
cients. Since the true deflection form represents a con- 
dition of minimum potential energy, any different form 
must represent a higher energy stage so that the 
buckling force must be greater than the true value. 





MATHEMATICAL ANALYSIS 


Consider a plate loaded in bending and compression at the plate middle plane. (See Fig. 1.) 


Boundary Conditions 


(1) w=Oonx=0,a; y=0,) 
(2) Vw=0O0onxs =0,a; y= 0,0 


The deflection surface is assumed to be of the form 


n= @ 


w = sin (rx/d) oa, sin (nry/b) (1) 


n=1 


where \ is the half wave length in the x direction. 


Each term of Eq. (1) satisfies boundary conditions 1 and 2. 


To determine the buckling force, the variation of the internal bending strain energy per half wave length, 6V, 
is equated to the variation of the external work. From reference 1 


terk (3 oy o*w ( ow \'} 
aie | Cel es - ———(——) |ode dy , 
; : J I (\ dx? * oy? ; "Lax? dy? dx oy) If” ( 








6ST = (1/DSL SN, (Ow/dy)? + N,(dw/dx)2] dx dy (3) 
Since 
w = sin (rx/d) > a, sin (nry/b) (1) 
n=l 
N, = No[1l — (2y/d)] 
Eqs. (2) and (3) may be evaluated in terms of a, 
5V = (Dbdw*/8) DO an?[(1/d?) + (2?/b?)}? (4) 
n=l 
(n = 2) odd 
ae ab n=©o i=@ F ; ™ 
67 = >: n'a,” + — No = = : (5) 
n=1 n=1 j=1 (n? =e 
Let 
, k Tr kyr? 
r/ b=8, - - = N ’ — = i Io 
b? . b? 
At the onset of buckling 6V = 67. Hence, 
co nm = 6 9 Noo wren A 
i a,(1 + Bn)? = kyBt DY a,’n® + ; — a Lb; —= = - (6) 
n=1 n=1 mw n=l i=l (nm? — c*)? 








(n + 1) odd. 
Ok,/ Oa, = 0 for k, a constant. 


[(1 + m?g*)? — k,Bin?]a, — 


(n + 1) odd. Solving for a,, 


i=o 


Eq. (6) is then minimized with respect to a, by differentiating with respect to a,, and setting 


= 0 ( 


168°k, > amt 


c tare = €) 


(3?/1687k,)a, = [Xo ayjin/(n* — 77)*]/[(11 + B?n?)? — k,Bin?] (7a) 


t=1 


(n + 1) odd. Eq. (7a) indicates that the odd Fourier 
coefficients may be given in terms of the even Fourier 
coefficients. Expressed in matrix notation, Eq. (7a) 
becomes 


1 (doaa) = [A }(Geven) l (x) 
yI (even) - [B (@oaa) f 


7 = unit matrix 


where () denotes a column, [] denotes a square matrix, 
and y equals the eigenvalue 7°/1687k,. Combining 
Eqs. (8) 

Ty*(Qoaa) =_ [A ][B ](doaa) = [C](Goaa) (9) 


where A is the square matrix of the even Fourier coeffi- 
cients, B of the odd coefficients, and C the matrix of the 
product of A and B. The advantage of this representa- 
tion becomes immediately evident, since Eq. (9) con- 
tains half the unknowns that would normally exist if 
both odd and even a,,s were unseparated. 


MATRIX ITERATION PROCEDURE 


If an exact solution is required, all of the a,’s must be 
retained. However, it is known that the higher modes 
contribute very little to the deflection form so that it is 
possible to reject some of them without altering the 
accuracy of the solution appreciably. To ascertain the 
number of coefficients required, the mode of buckling 
may be determined by an iteration method for finding 
the maximum eigenvalue and corresponding eigenvector 
of a matrix as outlined in reference 7. 

Initially, it is assumed that only four nonvanishing 
Fourier coefficients are necessary to determine k, 
within engineering accuracy. Setting 8 = land k, = 1, 
Eq. (7a) was expanded as shown below: 


ya; = (1/3)(0.222222a2 + 0.017777a4) 

yd, = (1/21)(0.222222a, + 0.240a;) | (10) 
ya3 = (1/91)(0.240a2 + 0.244898a4) 
yay = (1/273)(0.017777a, + 0.244898a3) 


where y = 77/1687. 


Operating on Eq. (11) in this manner 





Initial @) 
Matrix Guess x 10-* 
7,842,430 8,519,050 | 28 | 226,518 | 27.7118 
280,849 325.621 | 1 1 | 1 


In matrix form 
ay 7407.4 992.56 a2 
1057 = - , " 
. (“’) fosg 736-269. A (“:) 
., (a2 1058 .2 1142.85 | fa 
107 { “} = i ca 
’ (“) Iogpeotdl 89. fn (“’) 
Combining these two matrix equations 


(y108)27 (“’) _ i ty (“’) a) 

; a3 280,849 325,621 Az 

Eq. (11) is in the form required for matrix iteration. 
This method consists in taking an arbitrarily assumed 
set of values for a; and a3 and inserting them in the left- 
hand side of Eq. (11). A matrix multiplication will 
yield an improved column of a's. This procedure is 
repeated on the new column and is continued until the 
ratio of corresponding elements of the preceding column 
are not changed by further iteration. The maximum 
eigenvalue y,,,7 is then the ratio of the last column of 
a’s to the preceding columns. The a’s themselves are 
determined within a factor of proportionality. 

Svmbolically, this system may be represented as 


[A ](a,) = y(a,) 


where A is a known square matrix, a; is an unknown 
column, and y is the eigenvalue. Assume numerical 
values for the a’s. After a matrix multiplication, a 
new column of a's is obtained. 


[A ](a;”) = (a,) 


Operating on column (a;) in the same manner, a third 
column is obtained, etc. After 7 — 1 iterations 


[A |(a; “)) = (a;) 


J) (j-—1) 


( 
ay, (a; = 7 


where y is a constant; then y is the maximum eigen- 
value and the operation is completed. 





xe vill adn = 
226,313 | 27.8544 | 226,976 | 27.8546 
8124.95 | 1 =| 8,148.26 | 1 





L 


de 


Si 


ting 


(7a) 


ion, 
ned 
eft- 
will 
> iS 
the 
mn 
um 
| of 
are 


wn 
cal 


rd 


ELASTIC STABILITY OF FLAT RECTANGULAR PLATES 275 


v 


The second column in the first iteration is obtained by dividing through by 8,156.41. This process is called nor- 
malizing, and is performed to make the comparison of the iterated columns more convenient. Hence, 


Ymar2 = 8,148,260 X 10-", — a3/a, = 0.3590, — asa, = 0 
The same procedure was repeated using six nonvanishing coefficients. For this case, it was determined that 
Ymaz.. = 8,150,100 K 10-"°, a3/a, = 0.3586, a;/a, = 0.517 X 10-3 


A comparison of the two sets of results indicate that a; and ag may be neglected without decreasing the accuracy 
of the solution. Therefore, only a), a2, a3, and a4 were retained for the remainder of the analysis. Expanding Eq. 


(7a), 


ya, = (0.222222a. + 0.017777a4)/[(1 + 82)? — k,B*] 

ya = (0.222222a, + 0.240a3)/[(1 + 48%)? — 4k,B*] | (12) 
ya3 = (0.240a2 + 0.244898a4)/[(1 + 96%)? — 9k, B*] ae 
ya, = (0.017777a, + 0.244898a3)/[(1 + 168%)? — 16k, 8*] 


Let 
D = (1 + B)? — k, B64, E = (1 + 48")? — 4k,6' 
F = (1 + 987)? — 9k,6', G = (1 + 1687)? — 16k,6' 


and write Eqs. (12) in matrix form 


0.22222? 0.01777 —_—-(0.22222)(0.240)2 (0.01777) (0.244899) 
DE" DG =.lCU”.”C”C:*C~«U 
(0) ; () 
a3 (0.22222) (0.240) (0.01777) (0.244898) 0.240 0.244898 a3 
FE . FG FE + CD 


The stability determinant is obtained by subtracting y* from the main diagonal of the above matrix. If this 
determinant is set equal to zero, a polynominal in y is obtained of the form 














ae oe [oR (0.01777)? (0.240)? (0.244898)? ] 0.002399 (13) 
DE DG FE FG DEFG 
Since 
y= 1”/168°k, 
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Fic. 2. Instability coefficients for flat plates in combined bending and compression with simply supported edges. 








The largest root of the above biquadratic will yield the smallest value for &,. 


values of k, are plotted in Fig. (2). 


DISCUSSION 


The effect of a compression load on the instability 
wave form of a rectangular plate in bending is clearly 
shown by Fig. 2. Increasing the compression has the 
marked effect of increasing the wave length or decreas- 
ing the number of half waves in the direction of the 
bending stress. One interesting case is illustrated by 
the curve for k, = 1.25. For this condition the plate 
buckles in the direction of the bending stress with one 
half wave for a/b < 1.3, with two half waves for 1.3 < 
a/b < 2.4, and back to one half wave for 2.4 <a/b < 
2.9. 

The usual method for presenting the buckling of 
plates under combined loadings is through the use of 
interaction curves—that is, curves that determine the 
value of one stress required to produce buckling when 
a given value of another stress is also acting. Usually, 
the applied stresses are given as nondimensional stress 
ratios, R—i.e., a ratio of the stress acting to the stress 
required to produce buckling if that stress were acting 
alone. Interaction curves were constructed from Fig. 2 
and presented in Fig. 3. 

Interaction equations are usually assumed to be of 
the type 


RY? + Re + eee = ] 


where p and g are constant exponents chosen to fit 
experimental results. Further, this equation implies 
that any positive fraction of the critical stress of one 
type reduces the amount of bending of another type of 
stress required to produce instability. 

The present analysis indicates that for the problem 
investigated no single set of constant exponents (), q) 
exists which could represent the results correctly. It 
can be seen from Fig. 3 that the interaction curves are 
a function of the plate geometry, so that the exponents 
(p, g) cannot be independent of the plate aspect ratio 
as assumed in ANC-5. 

In addition, the interaction curves of Fig. 3 indicate 
that compressive interaction below the critical compres- 
sive stress has a small effect on the critical bending 
stress for long plates—i.e., for plates with length to 
width ratios greater than 2. However, compression is 
an important factor in the stability of short plates and 
should not be neglected in structural design. 

It is interesting to note that the matrix iteration 
method described in this paper was applied to the 
buckling of rectangular plates in shear and directed 
stress and failed to produce convergence. The reason 
for this difficulty was traced to the fact that the two 
smallest roots of the stability determinant were equal 
and opposite in sign. Physically, this means that the 
critical shear stress is independent of the direction of 
the applied shearing force. In order to obtain con- 
vergence, the geometric mean was taken of two succeed- 
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Fic. 3. Interaction curves for flat plates in bending and trans- 


verse compression at various aspect ratios. 


ing iterated columns to form a third column. After 
repeated iteration 0. this mean column convergence was 
attained. 

The same situation exists in the problem under in- 
vestigation: i.e., the buckling stress is independent of 
the direction of the bending moment so that the first 
two eigenvalues are equal but of opposite sign. To 
avoid the convergence difficulty described above, the 
matrix was reduced to a form where iteration produced 
the square of the eigenvalue. Thus convergence was 
assured, and, at the same time, the order of the matrix 
was reduced by a factor of 1/2. 
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Reduction Factors for Reinforced Monocoque 
Cylinders Subjected to Concentrated 
Radial Loads’ 


N. J. HOFF,* VITO L. SALERNO,t ann BRUNO A. BOLEY** 
Polytechnic Institute of Brooklyn 


SUMMARY 


Diagrams are presented from which the maximum bending 
moment can be calculated in a ring frame of a circular reinforced 
monocoque cylinder subjected to a concentrated radial load. A 
second set of diagrams is given for the determination of the maxi- 
mum shear stress occurring in the sheet covering of the cylinder 
in the vicinity of such a load. The curves are drawn on the basis 
of a theory that does not make use of the conventional simplifying 
assumptions of the engineering theory of bending. The results 
of this more accurate theory differ radically from those obtained 
from the conventional analysis. Reference is made to publica- 
tions describing experiments that corroborate the theory pre- 
sented. 


INTRODUCTION 


__gemaaagenaingy CARRIED OUT in the last 4 years have 
shown that the conventional theory of bending, 
when applied to the calculation of the bending moments 
in the ring frames and the shear stresses in the sheet 
covering of reinforced monocoque fuselages subjected to 
concentrated loads, gives values differing considerably 
from those measuredinexperiments. Maximum bending 
moments amounting to one-third of those obtained by 
conventional analysis and maximum shear stresses that 
are ten times those corresponding to engineering theory 
are common occurrences in the fuselages of modern 
transport planes. The reason for this discrepancy is the 
inaccuracy of the simplifying assumptions upon which 
the engineering theory of bending is based. In this 
theory it is assumed that the shape of the cross section 
of the beam remains unchanged when the loads are 
applied. Moreover, plane cross sections perpendicular 
to the longitudinal axis of the unloaded beam are 
assumed to remain plane and perpendicular to the 
axis after bending takes place. Centuries of engineering 
experience have shown these assumptions to lead to 
accurate and reliable results in the case of the solid cross 
sections customary in structural engineering. How- 
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ever, they are unsuitable as a basis for calculations in 
the case of the thin-walled and flexible reinforced 
monocoque fuselages of modern airplanes. When a con- 
centrated load is applied to a ring frame of such a fuse- 
lage, the ring frame is distorted both in its plane and 
out of its plane. These distortions result in a radical 
change of the stress distribution in the monocoque 
cylinder, and the magnitude of the changes depends 
upon the flexibility of the ring frame. Theoretical 
calculations showed that the stress distribution is in- 
fluenced primarily by a nondimensional structural 
parameter A defined as 


A = 6rt’/I,L? 


where r is the radius of the cylinder, ¢’ is the effective 
thickness of its material resisting normal stresses in the 
axial direction, J, is the moment of inertia of the cross 
section of the ring frame, and L is the distance between 
adjacent ring frames. The effective thickness ?’ can 
be calculated as the total cross-sectional area of all the 
stringers plus the effective sheet cross section, the whole 
quantity divided by the circumference of the cylinder. 
The parameter is essentially the ratio of the moment of 
inertia of the total fuselage cross section to the moment 
of inertia of a ring-frame cross section, multiplied by 
(r/L)*. A second and somewhat less important non- 
dimensional parameter is B, which can be written in the 
form 


B = 6(E/G)(t'/t)(r/L)? 


where E is Young’s modulus, G is the shear modulus, 
and ¢ is the thickness of the sheet covering. Physically, 
this parameter represents the ratio of the total cross- 
sectional area of the reinforced monocoque cylinder to 
the cross-sectional area of the sheet covering, the whole 
quantity multiplied by (r/L)?. 

Figs. 1 and 2 show typical shear stress distributions 
around the circumference of the fuselage for different 
values of the eccentricity parameter 7. This parameter 
is the ratio of the distance e from the skin median line 
to the centroid of the ring cross section and the radius 
r of the cylinder, as shown by the equation 


n=e/r 
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The shear stress distribution according to conventional 
engineering theory is also shown. It can be seen that 
the difference between the conventional and the more 
accurate theory is considerable and that the difference 
is greater for higher values of A. It may be mentioned 
that A = 10’ corresponds approximately to the fuselage 
of a modern transport plane, while A = 10* represents 
a fuselage of small diameter. B usually varies from 50 
to 400. Figs. 3 and 4 give similar information regarding 
the bending moment distribution in the loaded ring 
frame. Again, the difference between the conventional 
and the more rigorous theory is greater for the larger 
value of A. 

‘alues of two additional parameters—namely, C and 
D—are also noted on the diagrams. They characterize 
the extensional and shearing rigidities, respectively, of 
the ring frame. Their definitions follow: 


C = I,/A,r’, D = A,/A* 


a 


o 


In these equations A, is the cross-sectional area of the 
ring and A®* is its effective shear area. The latter is 
defined in such a manner that its product by Young's 
modulus is equal to the effective shearing rigidity of the 
ring. For thin-walled sections A* is approximately 
equal to (G/£) times the cross-sectional area of the web 
of the section. The values of C and D noted on the 
diagrams correspond to average design. 


Basic CHARTS 


The basic charts recommended for use in structural 
analysis are presented in Figs. 5 and 6. Engineering 
theory gives, for the maximum shear stress in the sheet 
cavering of a reinforced monocoque cylinder subjected 
to a concentrated load P, 


Tmax = P/(zrt) 


provided the monocoque is supported as a cantilever. 
The stress concentration caused by the deviation of the 
behavior of the cylinder from the assumptions under- 
lying the engineering theory can be taken into account 
through the use of a stress concentration factor u. The 
factor is defined as the number by which the maximum 
shear stress in the sheet covering according to the con- 
ventional calculation must be multiplied in order to 
obtain the maximum shear stress resulting from the 
more accurate analysis. The maximum shear stress is 
given, therefore, by the equation 


T maze = pP/(xrt) 


Calculated values of the shear stress concentration 
factor are plotted in Fig. 5 against the parameter A. 
The various curves correspond to different values of 
the parameters B, C, D, and». When the actual values 
of these parameters differ from those used in preparing 
the diagram, interpolation or extrapolation is facili- 
tated by the additional diagrams presented in a later 


section. 
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The location of the maximum value of the shear stress 
varies little. Approximate values of the angle at which 
the maximum occurs are listed in the diagram. 

Calculations showed that the maximum shear con- 
centration factor is obtained in a cantilever monocoque 
cylinder in the case when there is just one ring field in 
the cylinder. A ring field is defined as the portion of 
the monocoque extending between adjacent ring frames, 
The one ring-field cantilever monocoque therefore con- 
tains the loaded ring, the sheet covering, and the 
stringers only. It has no rings without loads. The ef- 
fect of the incorporation of additional unloaded rings 
between the loaded end ring and the rigidly fixed end of 
the cantilever is discussed in the next section. 

The maximum bending moment in the ring frame 
always occurs at the point where the concentrated load 
is applied. Engineering theory gives 


M maz. = (3/4m)Pr 


The actual behavior of the fuselage can be taken into 
account through the use of a moment reduction factor \, 
which is defined as the number by which the maximum 
bending moment according to conventional calcula- 
tions must be multiplied in order to obtain the maxi- 
mum bending moment of the more rigorous analysis. 
Consequently, the maximum bending moment is 


Maze = 3/4e)Pr 


Values of the moment reduction factor are given in 
Fig. 6 for different values of the parameters B, C, D, 
and 7. When the actual values of the parameters in a 
monocoque fuselage differ from those listed in the figure, 
interpolation or extrapolation is facilitated by the dia- 
grams presented in a later section. 

The values of the moment reduction factor shown in 
Fig. 6 were calculated for a monocoque cylinder that has 
an extremely large number of unloaded rings. The 
theory indicated that the maximum bending moment 
in the loaded ring is greater in such a case than in 
monocoques having a small number of unloaded rings. 
The analysis was therefore carried out for the semi- 
infinite cantilever monocoque, which extends from the 
loaded ring to infinity in one direction. The rigidly 
fixed end of the cantilever is at infinity. The effect of 
a finite number of ring fields is discussed in the next 


section. 


THE EFFECT OF THE ACTUAL LENGTH OF THE 
CANTILEVER 


The use of the diagrams of Figs. 5 and 6 always leads 
to safe design. When, however, a more accurate know'- 
edge of the values of the maximum shear stress and the 
maximum bending moment are required, use can be 
made of the information contained in Figs. 7 and 5. 
In Fig. 7 the variation of the shear concentration factor 
with the number of fields contained in the cantilever is 
shown. Curves are presented for a rather rigid struc- 
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ture (4 = 10*) and for a more flexible one (A = 107). 
The influence of different values of the parameters B, 
C, D, and 7» is also indicated. 

It can be seen that the shear concentration factor 
curve drops rather rapidly when the number of fields is 
increased, provided A is large. The variation in the 
shear stress concentration factor is insignificant when 
It appears that the value yu differs little 
the semi-infinite 


A is small. 
from its value corresponding to 
cylinder when the number of fields is greater than three. 

Fig. S is a similar plot of the values of the moment 
reduction factor. However, the curves reveal that the 
maximum bending moment increases when the number 
of fields contained in the cantilever is increased, while 
the preceding figure showed the maximum shear stress 


to decrease. 


THE EFFECT OF THE EXTENSIONAL AND SHEARING 
RIGIDITIES OF THE RINGS 


The stress distribution in a monocoque fuselage is in- 
fluenced not only by the bending rigidity of the ring 
frames but also by their extensional and shearing 
rigidities. The last two quantities are characterized by 
the parameters C and D. 

In order to facilitate the estimation of the change in 
the values of the shear stress concentration and moment 
reduction factors caused by variations in C and D, Figs. 
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9 to 12 are presented. Figs. 9 and 10 show the effects 
of the two parameters upon the shear stress concentra- 
Both diagrams were calculated for a 
10’). The variation 

that is, for small 


tion factor. 
relatively flexible cylinder (A = 
for relatively inflexible cylinders 
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values of A—-is much smaller. Similar statements can 
also be made in regard to the moment reduction factor 
whose variation with C and D is shown in Figs. 11 and 
12. 

In actual present-day design, C and D vary compara- 
tively little. The customary values are about 0.0002- 
0.0008 for C and 5-7 for D. In view of this slight varia- 
tion, it can be said that the effect of D upon the shear 
stress concentration and moment reduction factors is 
insignificant and that of C is generally small. 
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SHEAR CARRIED BY THE STRINGERS 


In all the preceding calculations it was assumed that 
the entire shear force was transmitted through the 
reinforced monocoque cylinder by means of shear 
stresses in the sheet covering. The stringers were 
assumed to be subjected to the normal stresses that 
equilibrate the moment of the load. Naturally, each 
stringer must have a finite bending and shearing rigidity, 
since otherwise it would buckle when subjected to 
compression. By virtue of their bending and shearing 
rigidities, the stringers are also capable of transmitting 
shear forces. Since any shear force taken by the string- 
ers relieves the sheet covering and, consequently, de- 
creases the shear stress concentration, a study was made 
of this effect, and the results obtained are presented 
here. 

In Fig. 13 the percentage of shear force taken by the 
stringer is plotted against the parameter B. (The dia- 
gram is valid for any value of A.) The type of loading 
shown in the figure is rather artificial, but it gives an 
indication of the effect sought. The two additional 
parameters H, and H, are defined in the following 
manner: 


H, = (mI aq.) /(2mt’r*) 
H, = (mI an) /(22t'r*) 


where /,,4. is the moment of inertia of a stringer for 
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bending in the radial direction, /,,,, is the moment of i io” ay ny — 
—— inertia for bending in the tangential direction, and m is A 
the number of stringers in the cross section. FIG.15 EFFECT OF STRINGER BENDING ON 
In Fig. 14 the same kind of information is presented ' SHEAR CONCENTRATION FACTOR 
for a somewhat different type of loading. The cylinder 
is now assumed to extend to infinity in both directions 
from the loaded ring. As it may be seen from a com- 
parison of the two diagrams, the percentage of the shear oy 
force taken by the stringers is not affected much by the bitin (") os 
difference in the loading conditions. 
As in actual present-day fuselages the values of the } | 2P 2P 
parameters H, and H, vary between 10~* and 10~, the —— — 
shear force carried by the stringers of actual fuselages } | —aauaieen 
I is small. Similarly, the effect of the bending of the 20| | = / 
stringers upon the moment reduction factor was found | |p op p op 


to be negligible. On the other hand, the shear stress 
concentration is reduced considerably by the ability | 
of the stringers to carry shear when the value of the a a 
parameter A is large. This can be seen from Fig. 15. | 











DIFFERENT CONDITIONS OF LOADING 





The information given enables the structural de- 
signer to calculate with little effort and with a reason- 
able degree of accuracy the bending moment in the 
loaded ring and the maximum shear stress in the sheet 
covering of a cantilever reinforced monocoque cylinder 
subjected to a concentrated radial load applied to the 
ring at the free end of the cantilever. Obviously, many 
other loading conditions occur in structural design for 
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FIG.16 SHEAR CONCENTRATION FACTOR FOR TWO FIELD 
SIMPLY SUPPORTED AND FIXED END MONOCOQUES 
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which information is also desired. Basically, a concen- 
trated load introduced into a cylinder is transmitted 
either in one direction or in two directions along the 
axis of the cylinder. For instance, in the fuselage of an 
airplane, a concentrated load may be equilibrated by the 
reactions of the fuselage at the wing attachments, or 
part of the load may be balanced by the aerodynamic 
and inertia forces acting upon the tail surfaces. When 
the former situation exists, the diagrams shown in this 
paper can be used with a reasonable degree of confi- 
dence. When the load is transmitted equally in two 
directions, the curves shown in Figs. 16 and 17 can be 
consulted. Fig. 16 shows the shear stress concentration 
factor corresponding to two-field cylinders whose ends 
are either simply supported or rigidly fixed. Compari- 
son of these two curves with those contained in Fig. 5 
reveals that the shear stress concentration decreases 
when the shear force is transmitted in two, rather than 
in one, direction. The two-field cylinder was chosen 
for this demonstration because it represents the most 
dangerous case of loading. The shear stress concentra- 
tion is smaller when the number of fields is greater. 

The moment reduction factors shown in Fig. 17 were 
calculated for the infinite cylinder. The infinite cylinder 
is defined as one that extends to infinity in both direc- 
tions from the loaded ring. Shorter cylinders give 
smaller moment reduction factors. It is noted that the 
simply supported and the rigidly fixed cases are not 
siown separately in this diagram. Both cases give the 
same moment reduction factors, since details of the 
distribution of the reactions at infinity cannot influence 
the stress distribution at the loaded ring. 

It can be observed that the moment reduction factor 
values in this case do not differ greatly from those 


presented in Fig. 6. 


HISTORY OF THE PROBLEM AND SOURCES OF 
ADDITIONAL INFORMATION 


The first published treatment of this problem was 
given by Wignot, Combs, and Ensrud,' of the Lockheed 
Aircraft Corporation, in May, 1944. They developed 
a differential equation for a somewhat simplified prob- 
lem and calculated shear stress and bending moment 
curves for many cases of loading and structural design. 
Their theory necessitated the use of an empirical factor 
and gave reasonable agreement with experiments car- 
ried out later. In June of the same year, Hoff? pre- 
sented a strain energy solution of the same problem 
using infinite trigonometric series for describing the 
shear flow distribution in the various ring fields. In 
the discussion of the paper, Newton* objected to the 
omission of the extensional and shearing strain energy 
terms for the rings. 

In the author's closure,‘ Hoff expressed his belief that 
the basic theory ought to be checked by experiment 
before it should be generalized to include the effects of 
the extensions and shearing of the ring frames. Such an 
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FIG. 17 MOMENT REDUCTION FACTOR FOR 
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experimental verification was obtained by the N.A.C.A. 
and published by Kuhn, Duberg, and Griffith’ in 
December, 1945. The N.A.C.A. tests were carried out 
with cylinders consisting of rings and sheet covering. 
The sheet was not reinforced by stringers. Strain-gage 
readings were in excellent agreement with the results of 
Hoff's theory. 

The extensional and shearing strain energy stored in 
the rings, the eccentricity of the ring attachment, and 
loading by concentrated moments and_ tangential 
forces were considered in a paper published by Beskin*® 
in June, 1946. In Beskin’s treatment of the problem, 
just as in reference 2, the calculations were carried out 
by the strain energy method with the aid of Fourier 
series. Beskin’s paper contains a large number of 
shear and bending moment curves. A different deriva- 
tion of expressions equivalent to those used by Hoff 
was given by Duberg and Kempner’ in 1947. The two 
authors considered radial and tangential loads, as well 
as concentrated moments, and prepared a large num- 
ber of shear stress and bending moment diagrams in a 
second publication.® 

In Great Britain, Goodey® presented a rather elabo- 
rate theoretical treatment of the problem in November, 
1946. <A different theoretical approach to the same 
problem can be found in a paper by Cicala'® published 
in September, 1946. Cicala considers each stringer 
separately, assumes the shear stress constant in each 
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panel of sheet, and develops and solves finite difference Auxiliary Properties (see body of paper for definitions) 


equations. The last two authors have not presented = 0.0234 
diagrams for practical applications. 0.00029 


yas 
Il 


The effect of the bending rigidity of the stringers was 5 
first treated by Levy.'' He assumed that only the t! 0.046 in. 
stringer at whose intersection with the ring the external A 6.69 X 10° 
load is applied has a finite bending rigidity. His theory B 256 
was corroborated by an experiment that is also de- H, 1.81 & 10-5 
scribed in the paper. Another solution of the same H, = 1.16 X 10-5 
problem was given by Hoff'* in 1947. He distributed 
the bending rigidity of the stringers uniformly around 


Consider a load P = 1,000 Ibs. at the end of a three- 
field cantilever. 

In Fig. 5, interpolation between Cases II and IV gives 
for the shear concentration factor n = 10.8 in the one- 


the circumference. 
Additional investigations were undertaken at the 
Polytechnic Institute of Brooklyn under a research 


contract with the office of Naval Research of the United field case. Hence, 
States Navy. Most of the diagrams presented in this r maz. (conventional= 1,000/zrt = 165 Ibs. per sq .in. 
paper are the results of those investigations. Details , (present theory) = 10.8 X 165 = 
c . . . mar i — «< e = 
of the theory and a comparison with experiments were 1,780 Ibs. per sq.in. 


published in January, 1949."% 
The purpose of the experimental part of the Poly- In Fig. 6, interpolation between the curves corre- 
technic Institute of Brooklyn investigations was to ex- sponding to B = 250, C = 0.0005, D = 5,» = 0, and 
tend the N.A.C.A. results in several directions. One B = 250, C = 0.0005, D = 5, » = 0.075 yields, 
extension was in the values of the parameters A and B-_ for the semi-infinite cantilever, ) = 0.485. Conse- 
covered. While in the N.A.C.A. tests A varied from quently, 
L815 to 155,000 and B from 4.0 to 18, in the Polytechnic - M mex. (conventional) = 3Pr/4x = 14,300 in.Ibs. 
Institute of Brooklyn test the maximum value of A Vv (present theory) = 0.485 X 14,300 = 6,930 
was 10° and that of B was 451. Tests were also carried ee ’ ; = 
out with rings attached to the cylinder externally and 
internally in order to check the effects of the eccentric- _ therefore, 
ity. Moreover, cylinders were constructed with and 
without stringers and with various values of the param- qe . 
eter /7, and H, in order to investigate the effects of the ne 65,700 Ibs. per sq.in. 
finite bending rigidity of the stringers. The agreement Omaz- (present theory) = 0.485 X oer re ; 
between theory and experiment was satisfactory. $1,800 the. per sq.tn. 


The maximum bending stress in the loaded ring is, 


O mare (Conventional) = 14,300/0.218 = 


The values calculated can be considered as safe 
NUMERICAL EXAMPLE approximations. Correction for the actual number of 
ring fields gives in the case of the shear stress (see the 


The fuselage to be investigated has the following 
full line corresponding to A = 10? in Fig. 7). 


geometric properties: 


Monocoque Cylinder hw = (9.5/13)10.8 = 7.89 


t (thickness of sheet) = 0.032 in. a nee eee " 305 It i 

r (radius of monocoque) = 60 in. +300 Ibs. per sq.in. 

L (ringframe spacing) = 18 in. In the case of the maximum moment the correction 
Stringers can be made with the aid of the full line of Fig. 8 (cor- 


: : responding to A = 10’): 
A x, (area of a stringer) = 0.133 in.? — aes ) 


I, (moment of inertia in radial direction) = 0.0282 X = (0.37/0.418) (0.485) = 0.428 
in.4 Omar. (present theory) = 0.428 (65,700) = 
/,(moment of inertia in tangential direction) = 28,100 Ibs. per sq.in. 


0.0181 in.4 


m (number of stringers) = 40 It can be seen from Figs. 9 and 11 that the effect of C 


‘ (the true value is C = 0.00029, while the curves used 
Ring frame were calculated for C = 0.0005) is insignificant. Fin- 
ally, Fig. 15 indicates the possibility of a further reduc- 
tion of the maximum shear stress in the skin by 
Z (section modulus of ring frame) = 0.218 in.* approximately 5 to 10 per cent because of the finite 
A,, (area of web of ring) = 0.17 in.” bending rigidity of the stringers (WH, = 1.81 X 10~° 
A* (effective shear area of ring) = 0.066 and H, = 1.16 X 107°). 


A, (area of one ring) = 0.330 in.? 
! (moment of inertia in its plane) = 0.350 in. 
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The Theory of Beam Columns with 
Continuous Support Applied to the Problem 
of Stringer Roll 


G. C. BEST* 
Consolidated Vultee Aircraft Corporation 


ABSTRACT 


In this paper the theory of beam columns of constant section 
with continuous support is developed, the mathematical pro- 
cedures following those already laid down by Timoshenko and 
other writers.') *} 8 All deflections are considered to be .within 
the proportional limit. Only the case of a moment at one end is 
considered, the principle of superposition then permitting its 
extension to the case of a moment at each end. Buckling loads 
are obtained from the criterion that moments become infinite. 
Results are applied to the problem of stringer roll by considering 
the vertical portion of the stringer as supplying continuous sup- 
port to the flange, which is considered as a beam column. An 
illustrative example is worked out. 


NOTATION 


P = axial load in elastically supported beam 
L = length of span between pin supports 

£ =a dimensionless coordinate = x/L 

7 = adimensionless coordinate = y/L 

EI = flexural rigidity of beam column 

k = spring constant of elastic support per inch length of beam 
I 


?.. = buckling load of continuously supported column of in- 
finite length = 2\/REI 

Lo. = wave length into which a continuously supported column 
of infinite length buckles under the load Po 

P, = buckling load of pin-ended column with no elastic 
support = r?EI/L? 

u =a dimensionless parameter = (1/4/2)L\/P/EI = 
(1/V/2)(L/j) = (#/V2)V P/Pe 

#® = adimensionless parameter = P./P 


a dimensionless parameter = L/0.5L 2. 


> 
Ul 


Basic THEORY 


I THE FOLLOWING ANALYSIS, deflections are considered 
to be sufficiently small to permit the substitution of 
d*y dx? for M/EI. Plasticity effects are neglected, and 
only beams of constant section are considered. The 
elastic support and also the deflections are considered 
to exist only in the plane of the applied moments. 
Details of the mathematical derivations are given in 
the Appendix. 

To obtain the equation of the beam column with con- 
tinuous support, it is only necessary to add the effect 
of the continuous support to the already existing equa- 
tions for the beam column. Measuring coordinates as 
indicated in Fig. 1, the equation of the beam column 
with transverse loading w(x) and end moments M, and 
M2 is: 


Received September 7, 1948. 
* Stress Analyst. 


M = M, — Rix + Lf w(x) dx dx — Py (1) 
0 Jo 


The moment due to the continuous support is equiv- 
alent to that caused by a loading —ky, where k is the 
spring constant of the support per inch length; hence, 
the equation of the beam column with continuous sup- 
port is: 


M= M — Rix - f I ky dx dx — Py (2) 
0 0 


Differentiating twice gives: 


@d?M/dx? = —ky — P(d*y/dx?*) (3) 
which, since d?y/dx? = \//EI, becomes 
(d?M /dx*) + (PM/EI) + ky = 0 (4) 


equivalent to reference 3. 
To shift to dimensionless coefficients, let & = x/L 
and » = v/L. Then 


d?M/dx*? = (d?M/dé*)(1/L*) 
and Eq. (4) becomes: 
(d?M/dé&*)(1/L?) + (P/EI)M + kLyn = 0 (5a) 
This can be written 
—kL'n = 2aM + (d*M/dé*) (5b) 


where a = PL?/2EI. This form will be useful later in 
determining deflections. Differentiating Eq. (5a) 
twice with respect to £ and noting that 


d?y/dx? = (dn/d#)(1/L) = M/EI 


gives: 


a1M 1 P\/d?M ,f M be 
(SG) + alas) seas ( EI ) =" 


(6) 


a‘M +4 (\(S) + (=) M=0 7) 
dé* EI /]\ dé EI} “ 
all coefficients now being dimensionless. Letting 5? = 


kL‘/EI gives 
(d‘M/dé&*) + 2a(d?M/dt?)b°?M = 0 (8) 





or 





or in symbolic notation 
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(D4 + 2aD? + b2)M = 0 (9) 


It should be noted at this point that for any fixed axial 
load P the coefficients a and } are constants independent 
of M. Consequently, Eq. (9) is a homogeneous linear 
differential equation and solutions may be superim- 
posed. Now let D = D®, transforming Eq. (9) into 


(D? + 2aD + B*)M = 0 (10) 
or 
D = -¢a + Va? — 3b? (11) 


There are three cases to consider 6? > a?, 6? = a*, and 
b? < a®. For 6? = a? or RL*/EI = (PL?/2EI)?, one 
obtains 


P, = 2VREI (12a) 


This corresponds to the critical buckling load* !' of a 
continuously supported column of infinite length 
buckling in waves of length® !° 


L. = 24 VEI/k (12b) 
It is convenient to let 
w = P,/P = 2VREI/P (12c) 
then 6 = aw. Eq. (11) then becomes: 
D= —a + aV1 — w? (13) 


The three cases requiring consideration become: I, 
» >1; I,w = 1; HI,w <1. 


Case I (w > 1) 


, / , 7? 
In this case V1 — w? is complex; hence, Eq. (13) 
becomes: 


D = -a +aVwt—1 i = —re*¥ (14a) 
where 
i= V-1 
r= aw 
P _ 
g = tan“! + Vow’ — 1 
also let 
u = (1/V2)(L/j) (14b) 
where j = WV El/P following the notation of reference 


lsothata = wu? andr = wu?. The roots of Eq. (9) are, 
therefore, since D? = 0: 


MAY, 1949 


+ isin .) 


6 
6 
I 


7 
~ 


f/f . 9 i / 
gy rev Ol" = iV r (cos 


4 4 ° gy 
+Y “(sin > = 1 cos .) = +(a@ + §8) 
where 
wou 7) pr 1 =~ i/e) w— | 
a= Vary 3 “e@ . 
D r 1 + (1/w) lo +1 
¢ - _ 
B= %N r C08; = VV wun, > = ¥ \ 9 


(l4c) 
The general solution of Eq. (9) is, therefore, 


M = e“(A sin BE + Bcos BE) + 
e~*(C sin BE + D cos BE) (15) 


For the case 14, = 0 and M, = 1 one obtains A = C 
and B = —D, details being given in the Appendix. 
Hence, Eq. (15) may be written: 

M = 2A sin Bé cosh at + 2B cos BE sinh aé (16) 


End conditions [Eqs. (25)] determine the constants A 
and B. The following results are obtained: 


A = ko ead ek, 
B= ky + eke 


where — 
€ = 1/Vw? — 1, ky =d/(e? + a2), be = e/(e2? + a?) 


and 


d = 2 cos B sinh a 
e = 2sin B cosha 


The point at which maximum moments occur may be 
determined from the equation: 


tan B§ = K/tanh aé (17) 
where K = (aB + BA)/(8B — aA), which may be 
solved by iteration. 

Case II (w = 1) 


In this case both roots of Eq. (10) are equal. Eq. (14) 


yields D = —a; hence, the roots of Eq. (9), which are 


both double, are + i~/a = +iu. 
The general solution is, therefore, 


M=Asinui + Boos ué + 
&(C sin ué + Decos ué) (18) 


As shown in the Appendix, end conditions again deter- 


mine the constants for the case 1, = O and Mz = 1. 
Results are: 


A = (2sin u — ucos u)/2 sin? u 
B=zCz2=0 


D = u/2 sin u 


giving: 








= 18) 


l4c) 


(15) 


dix. 


16) 


be 


be 


ire 
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M =A sin ué + Dé cos ué (19) 


it being possible to determine the position of the maxi- 
mum moment from 


tan ué = G/ué (20) 


where G = (Au + D)/D, solving by iteration. 


Case III (w < 1) 

For this case Eq. (10) has two negative roots, since 
V1l—w?< 1. Let these be —8,?, —6.”. Then, from 
Eq. ( 13), 

—B, = —a(l — V1— w?) 
— B.? —a(1 + V1 — wo?) 


or 
6 =uVi1— V1 — w? 
Jf — — 
Be = uVi+ V1 — wo 


The roots of Eq. (9) are, therefore, +781, +762, the 
general solution becoming: 


M=A sin Bit + Boos Byé + 
C sin Bt + D cos Bt (21) 


For the case 1/,; = 0, 2 = 1, the constants in Eq. (21) 
become (see the Appendix): 


A = 1/(1 — y) sin B; 
B=D=2=0 
C = —y/(1 — y) sin Be 


where 
y = B2/B2 = (1+ V1 — w)/(1 — V1 — @) 
hence, 
M =A sin BE + Csin Bot (22) 


The position of the maximum moment may be deter- 
mined from 


HT cos BE = cos Boe 
where 


H = (sin B2/sin B;)/y(B2/:). (23) 


SUMMARY OF CASES 


Though the form of the analogous Eqs. (16), (19)» 
and (22) in Cases I, II, and III is different, Eqs. (16) 
and (22) nevertheless approach Eq. (19) as a limit as 
# = |,asisshown in the Appendix. The determination 
of the maximum moment can be accomplished, if exact 
answers are sought, by means of the corresponding 
Eqs. (17), (20), and (23), but the most practical pro- 
cedure is simply to plot \/ against £, this being advisable 
anyway to obtain an idea of the shape of the deflected 
beam. Actual deflections, if desired, can be obtained 
by inserting the final equation for moments into Eq. 
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(5b), and, of course, slopes may be obtained by differ- 
entiating this result. 

As FP increases, it is shown in the Appendix that as 
long as k is also increased so as to maintain Case I 
relationships, 1/ does not become infinite. For Case 
II, however, it is evident from Eq. (19) that for u = 
ma where m is a positive integer, the sin « = 0 and 
moments tend to infinity. This corresponds to the 
buckling load. Similarly, in Case III for either 8; or 
82 = mm, denominators approach zero and moments 
become infinite as the buckling load is approached. It 
is advisable, therefore, to compute first the buckling 
load before applying the above formulas, since results 
given by the above equations for / are meaningless for 
loads exceeding the critical. The buckling load may be 
computed from reference 8, which is also shown in this 
paper as Eq. (42). 


APPLICATION TO STRINGER ROLL 


By considering the stringer as theoretically composed 
of separate parts, see Fig. 2, where the shaded portion 
constitutes a cross section of the “‘column”’ involved and 
the vertical web and skin provide the ‘“‘continuous 
elastic support,’ it is possible to apply the results of 
the above theory to the problem of stringer roll. The 
determination of how much of the vertieal element to 
include with the flange of the ‘‘beam column”’ and the 
selection of points from which to estimate deflections 
of the ‘‘continuous elastic support”’ requires the exercise 
of judgment, especially since in the neighborhood of 
buckling loads slight variations in these measurements 
will cause large variations in final maximum moments. 
Reference should be made to test data, of course, as 
far as possible. When no end moments exist and the 
buckling load alone is sought, the above method is still 
applicable, but it is preferable in such cases to employ 
the more exact procedure of reference 7. 


EXAMPLE 


To illustrate the application of the method the fol- 
lowing problem is worked out. A stringer of section as 
shown in Fig. 2 is mounted so that the unsymmetrical 
flange is free between bulkheads, which are 22 in. apart. 
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The shaded portion is considered to be the beam, and 
the vertical web and skin form a system providing con- 
tinuous elastic restraint. Adjacent stringers are located 
at A and C, Fig. 2, the skin being considered fixed at 
these points and pinned at B. Consider a section | in. 
deep. The moment of inertia of the skin or web is 
1(0.125)*/12 = 0.0001625. E is asumed 10.5 X 10°. 


The deflection caused by a 1-Ib. force at D is: 


1 X (3.3)8 : (3.3)? x 
3X 10.5 X 10® & 0.0001625 2 
10 


———— - = 0.0152 
4 X 10.5 X 108 X 0.0001625 





Hence, k = 1/0.0152 = The J of the shaded 
portion is 0.0915 in.‘, and the area = 0.4 sq.in. 


65.7. 





P.. = 2VkEI = 
2V65.7 X 10.5 X 10° X 0.0915 = 15,800 Ibs. (12a) 
P., = 22,700 Ibs., m= 1 (42)8 


Assuming a compressive stress of 40,000 exists, P = 
0.4 X 40,000 = 16,000 Ibs. 


a | _ 16,000 on 
u = a 4 = —— oS Ee a os = ae 
V2) ~ 7/2 N10.5 X 10° X 0.0915 
P., 15,800 7 
wy = — = —— = 0,985 (Case III) 
P 16,000 
1 — w? = 0.029775, W1 — w? = 0.172555 


Vi — V1 — w = V0.827445 = 0.90970 


Vi4 V1 — 0 = V1.172555 = 1.08284 


— = 1.41708 


1.172555 
5 





re 1—-Vi— wt 0.82744 
1 Pens or ” 
= —23976 —- = 4+3.3976 

l-y¥ ey 


8; = 2 X 0.90970 = 1.8194, 
B. = 2 X 1.08284 = 2.1657 
sin 8, = 0.96927, sin 6B. = 0.82820 
A = —2.3976/0.96927 = —2.4736 
C = 3.3976/0.82820 = +4.1024 
Eq. (22) becomes 


M = —2.4736 sin 1.8194 — + 4.1024 sin 2.1657£ 


Results for various values of £ are given in Table 1. 
The maximum value of M for 1 in.lb. at end & = 1 is 
1.758. The resultant load in the shaded portion shifts 
through a distance of 0.5150 in. at the point of attach- 
ment of the jay stringer to the adjacent symmetrical I 
stringer. This causes an eccentric moment of 16,000 X 


AERONAUTICAL 


SCIENCES—MAY, 1949 
TABLE 1 : 
é Big Boe sin f; & sin Bot M 
0 0 
0.1 0.18194 0.21657 6.1809 0.2149 0.4340 
0.3 0.54582 0.64971 0.5191 0.6050 = 1.1977 
0.5 0.90970 1.08285 0.7893 0.8833 1.6711 
0.651 1.18463 1.41011 0.9264 0.9871 1.7580 
0.8 1.45552 1.73256 0.9934 0.9870 1.5917 
1.0000 


1.0 


0.5150 = 8,230 in.lbs. Three-fourths of this moment 
is assumed carried by other structure and one-fourth 
or 2,060 in.lbs. by the beam. Hence, the actual maxi- 
mum moment in the beam is 2,060 X 1.758 = 3,620 
in.lbs. giving a total stress of 

3,620 X 0.5775 


: = 62,800 Ibs. per sq.in. 
0.0915 


frot. = 40,000 + 


A ppendix 
More complete details of the mathematical proced- 
ures in the above are as follows: 
Case I 
Differentiating Eq. (15) twice and substituting in 
Eq. (5b) gives: 


—kL'n = Ae®™[(2a + a? — B*)sin BE + 2aB cos 
Bt] + Be“ [(2a + a? — 8?) cos BE — 2aB sin 
Bt] + Ce-“[(2a + a? — B*) sin BE — 2a8 cos 


BE] + De-™ x 


[(2a + a? — B*) cos BE + 2aB sin BE] (24) 


The énd conditions from which the constants may be 
determined are: 


For — = 0, M = M,and yn = 0\ 25) 
For — = 1, M = Meand n = Of - 
Eqs. (15) and (24) then yield, letting « = (2a + 
a? — B)/2aB = 1/Vw? — 1, 
M, == B + D 
O=A+e«eB-—-C+e&D 
M, = e*A sin B + e“Bcos B + e~*Csin 
8+e-Dcos Bs } (26) 


O = e*[esin 8 + cos BJA + e*[ecos 8 — sin B] | 
B+ e-*[esin B — cos BJC + 
e~“*[ecos 8 + sin BJD | 


This simplifies to: 


M,=B+D 
O=A-—-C+eM, 


M, = A(e* + e-*) sin B + B(e* — e-*) cosB+ (97) 
Mye-*(e sin B + cos 8)( 
O = eM, + eA cos B — e*B sin B — e-* X 
(A + €M;) cos B + e~*(.M, — B) sin B 
Letting 


m, = M, — Mye-*(e sin B + cos B) 
eM, + Mye-*(sin 8 — € cos B) 


mM, = 





340 
77 
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reduces Eq. (27) to 


II 


mM, eA + = (28) 


—m, = dA — eB’ 

which, solving for A and B, yields 
A = mkz — moki, B= km + kom. (29) 
Setting M@, = 0 and M2 = 1 gives Eq. (16). Differ- 


entiating Eq. (16) with respect to £ and equating to zero 
yields Eq. (17). 


Case II 
Differentiating Eq. (18) twice and inserting in Eq. 
(5b) gives: 


—kL*n = aA sin ut + aB cos ué + 
[aé sin ué + 2u cos ué]C + 
[aé cos ué — Qu sin ué]D (30) 


Consideration of end conditions (25) together with 
Eqs. (18) and (30) yields: 


M,=8B 

O= (wW/)B+C 
Me =Asinu+ Bcosu+ Csinu + Deosu (31) 
0=Asinu+ Bceosu+ 


[sin uw + (2/u) cos u]JC + 
[cos u — (2 sin u;u)|D 


Setting M, = 0 and M, = 1, this simplifies to Eq. 
(19), Eq. (20) being determined from Eq. (19) by dif- 
ferentiating and equating to zero as before. 


Case III 
Differentiating Eq. (21) twice and inserting in Eq: 
(5b) gives: 
—kL*n = (2a — B,")(A sin BE + B cos BE) + 
(2a — B.”)C(sin Bet + D cos Bet) (32) 


End conditions (25), together with Eqs. (21) and (32), 





letting 
2a— 8B? Be 1+ V1—a? on 
SS = _= 3: 
2a — B2* Bi" i x SE ae a 
gives 
M=B+D 
0=7B+D 
M, = A sin 8; + B cos 6, + C sin 2 + (34) 
D cos Bo : 
0 = (A sin 8; + B cos 8:1) + C sin B. + 
D cos Bo 
from which are obtained: 
B= M,/(1 — vy) 
D= (—y/1 —_ 7) My, (35) 
A = (M2, — Mi cos 6;)/(1 — y) sin pi 
C = [-y/(1 — y)][(M2 — Mi cos B)/sin Be] 


For M, = Oand M; = 1, Eq. (22) results, from which 
Eq. (23) may be obtained by differentiating and equat- 
ing to zero as before. 


LimITING Case (I = IT) 


That Case I has Case II as a limiting case as w = | 
with w > 1 can be shown as follows. Let w = 2h? + 1 
so thath = Oasw = 1. Then 


Pe 
a = uh, B=uVih?+1= 
d = 2 cos B sinh a = 2uh cos u 
e = 2sin Bcosha = 2sinu 


Vo? — 1) = 1/2h 


Therefore in the limit, when J1/; = 0 and MJ, = 1, 


e— ed 2sinu — 2uh cos u 2h) 
e? + da (2 sin u)? % 
2 sin u — uCcOs u 
(2 sin u)? 
_e+d [2 (sin u)/2h] + 2uh cos u 
et + d? - (2 sin u)? 7 


] 


2h (2 sin u) 





hence, from Eq. (16), 


stl, sin B§ cosh aé + 
of = = , cos Bé sinh at 
2h 2 sin u . 


? 
sin “ — u Cos u ug cos ué 
M= (7= —- ‘) sin wé + 


2 sin? u 2 sin u 


: 2 sin u — u cos u 
Me 
(2 sin u)? 


which is identical with Eq. (19) of Case IT. 


LimITING Case (III = IT) 


To show Case IT to be the limiting case for Case IIT 


= 1—/h? or V1 — w? = 


asw = 1 withw < 1, let w? 
Then 


8,2 = a(1 — h) 
8, = uV1 —h =ul[l—(1/2)h...J =u 
Bo? = a(1 + A) 
=uVl+h=u[l+ (1/2)h...J =u (36) 
sin Bi = sin[tu — (1/2)hiu] = 
sin &u cos (1/2)htu — sin(1/2)htu cos tu = 
sin &u — (1/2)héu cos Eu |} 


with similar formulas for sin ft, sin 6), and sin Be. 
Furthermore, 


ae 2 1 


sin u — (1 2)uh c cos u 


1 cos u a 
ae 
sin u sin u 


sin 2; 
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Inserting values for A and C in Eq. (22) gives: 
l sin Big oY sin Be Pi 


M = —_ 
~ (1—y) sin B; (1 — ¥) sin Bs 


l l l cos u 
ee aoe 1 + — uh | — x 
(1 — y) \sin w 2 sin u 
l 
(sin tu — —-héu cos tu) _ (2 Ysa - x 
2 1 — y/\sin u 
l cos u 
~~. ar sin &u + - su cos £u 
2 sin u 


Expanding and observing that y = #:°/f:> = 
(1 + 4)/(1 — h), so that (1 y)/QA — ¥) = -I/h, 
gives Eq. (19) of Case II. 

oo ) 


LIMITING CASE (w = 


If k and P both increase so that w > 1, the formulas 


of Case Lapply. For large values of w, a = 8 = u V w/2 


e* cos 8, e = e* sin 8, since for large values of a 
one obtains 2 sinh a ~ 2 cosh a = e*. Also € = 0. 


A = e-“sin Band B = e~* cos 8, and Eq. (16) 


andd = 


Hence, 
goes over into 


M = e~*'—® sin B sin BE + 
e~*"'—® cos B cos BE = 
e-atl—8) cm- 2 2 
cos B(1 §&) $1 (39) 


so that 7 does not approach © with P in Case I. 
If k is held fixed and P = 0, w = 


it is convenient to write 


To study 
Then 


© again. 


this case p = l/w. 


pb = a, and Eq. (11) becomes: 
D) = —pb + V p°b? — bh = —bd(p = Vp? = 
(40) 
For p = 0 the roots of Eq. (10) are +bi = a he(*/?)* 


as / 

and the roots of Eq. (9) become: +V 6 et =Vb 
‘ . / J is rege A 
ie Pt eV b ef er or =(V h os b 21), which 


yields the general solution: 


l l 
M = ev va A sin \ <t + Bcos \ 4 a 


— Db b 
e C sin V5 + D cos V5 (41) 


which corresponds, except for notation, with reference 
6—i.e., the case of continuous support but no axial 


load. 
LIMITING CASE (w = 0) 


If P is held constant and k = 0 then w = O and the 
Eqs. of Case III apply. Then §,;? = 0, B.? = 2a = 
(L/7)*, oF hts = L/j in the notation of reference 1. 
Also y = », and Eqs. (35) become B = 0, D = +M,, 
A = 0, C = [Mz — M, cos (L/j)]/sin (L/j), and Eq. 
(21) Adee identical with reference 2. 


RONAUTICAL 


SCIENCES MAY, 1949 


BUCKLING Loaps (CASE II) 


Since sin “ appears in the denominators of the coef- 
ficients A and D of Eq. (19) of Case II, the moments 
as u« = mm, where m is an integer. 
V2) L V P/EL for u and clearing radi- 
cals gives P = 2m°P,. Since in Case II P = P,, this 
becomes m? = P,,/2P, = \*. Hence, if \ is an integer, 
buckling occurs in Case II at the load P,. It can be 
shown®~!! that m (or \) is the number of half waves in 
the buckled form of the supported beam—i.e., for buck- 
ling at P., in Case II, Z must be an integral multiple 
of 0.5Z,, a result that might have been expected 
intuitively. 


tend to infinity 
Substituting (1 


BucKLING Loans (CAsE III) 


Similarly, in Case III it can be seen from Eq. (22) 


that moments again approach infinity as either sin A, 
— a oe / / ner 
or sin §. = 0. Setting w¥ 1+ sV 1 — w? = ma where 
+1 or —1 and m is an integer, substituting for u 
and w and clearing radicals gives 


P., = P,[m? + (A‘t/m?*)] (42) 


which is identical, except for notation, with reference 8. 
This equation may be used to determine the buckling 
load. In this formula m is an integer so chosen as to 
make P., a minimum, or m may be taken as the greatest 
integer in 


(1 + V1 + 4a2)/2 (43) 


obtained from reference 9. In obtaining Eq. (42) above 
for s either +1 or —1, only one such development is 


valid. For, in the process, the equation 


sV1 — w? = 2m°(P,/P) — 1 (44) 


is obtained, which is true for only one of the values s = 
+1 or s = —1 depending on whether 2m*P,/P is 
greater or less than unity. If 2m°?P,/P = 1, Case II is 
implied. However, if \* = m(m + 1), the beam may 
buckle into either m or m + 1 half waves,’ in which 
case, if 6; = m7, B2 will equal (m + 1), as can be 
shown by combining the relations of references 8 and 9. 

Holding k, EJ, and L constant and increasing P yields 
the result that 8: increases and (; decreases with increas- 
ing Pin Case III. This can be shown by substituting 
for « and w in the formulas for 6; or 62 and clearing 
radicals, yielding 


| ‘PF i ae 
P= “<a + ll P ( ) (49) 
Tr 4 P, B? 
where 7 = 1,2. Differentiating P with respect to 0; 


yields 
dP w? 
dB, (1+ V1 — w)? 


the + or — sign going withz = 2 or 1, respectively, the 
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bracketed expression accordingly becoming positive or 
negative. 
STIFFNESS AND CARRYOVER FACTORS 


To obtain stiffness and carryover factors it is con- 
venient to use a different sign convention. Positive 
rotation will be defined as counterclockwise and a posi- 
tive moment in a beam as that which tends to rotate a 
joint counterclockwise. End rotations may be obtained 
from Eqs. (24), (30), and (32). Let 


aga = (dy dx). =(90 = (dy dé)» = @ 
a, = (dn/dé); — 1 


then carryover factor 
C = —ao/ay (51) 
and stiffness correction factor, other end fixed, 

o = (L/4ET) [a1/(a1? — a?)] (52) 
where VM = aS and S = (4E//L)c. These formulas 
yield the following results: 

Case I 
C = 1/(¢ sinh a sin B — cosh a cos B) (53) 
where ¢ = (ak; + Bke)/(ake — Bh). 
o = [(u?/a? — 1)/8(ak, — Bki)][C/(1 — C?)] (54) 
Case II 
C = (sin u — ucosu)/[u — (sin u)(cos u)] (55) 


o = (u/2)({u — (sin u)(cosu)]/(u2 — sin? u) (56) 


Case III 


C = —(6, sin B, — BP sin 6) + 
(81 cos 81 sin Be — Be cos B. sin By) (57) 


V1 i =| B, cos 8; sin Be — Be cos Be sin p; | 
_—— = ~ Lw(1 — cos 6; cos B2) — sin B; sin Be 


(58) 


OTHER NOTATIONS 


The constants u, a, B, Bi, Bo, etc., in the foregoing 
can also be expressed conveniently as functions of p 
and A, where p = 1/w and A = L/0.5L., both being 
ratios of similar significant elements, therefore making 
excellent dimensionless parameters. In this notation 
one obtains 


a’? P.. 


t= PLP, wut = >Sent (47) 


u = rAVp (48) 


(49) 


B = rAVp + Vp? — 1 (50) 
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Loads on a Supersonic Wing Striking a 
Sharp-Edged Gust 


M. A. BIOT* 


Cornell Aeronautical Laboratory 


ABSTRACT 


This is a calculation of the chordwise lift distribution, total lift, 
and moment on a two-dimensional wing striking a sharp-edged 
gust at supersonic speed. A direct solution is established by 
considering a distribution of sources in a fluid at rest. An alter- 
nate method using Busemann’s conical flow is also shown to be 
applicable. The time history of the total lift and mid-chord 
moment is discussed. It is shown that the total lift increases 
with time and reaches a maximum that corresponds to the steady- 
state phase of the flow. The mid-chord moment goes through a 
maximum independent of the Mach Number if the latter value is 
larger than 4/2, while this maximum can become infinite for a 
range of Mach Numbers between 4/7 and 1 


(1) INTRODUCTION 


_ Pipes ATTENTION has been given lately to 
nonstationary flow problems of wings flying at 
supersonic speeds. Most of the work, however, has 
been concerned with the aerodynamic forces on an 
oscillating airfoil from the standpoint of flutter analysis. 
The problem of the wing hitting a sharp-edged gust is of 
a different nature and turns out to be actually much 
simpler than the oscillating airfoil problems. 

It is shown in section 2 that it may be treated by a 
distribution of sources of a simple type along the chord 
and that the pressure distribution may be derived by 
elementary methods. The procedure does not intro- 
duce a moving fluid but considers a fluid at rest in 
which nonstationary sources are distributed in a layer of 
variable extent. This point of view, which is closer to 
acoustics than to aerodynamics, is somewhat novel and 
seems to present advantages of simplicity and closeness 
to physical reality in certain categories of problems. 
The pressure distribution derived by this method is 
applied to the calculation of the time history of lift and 
moment on the wing in section 3. Particular attention 
is given to the value of the mid-chord moment, which 
starts from zero, rises to a maximum, and goes back to 
zero. The value of this maximum and related data is 
evaluated in section 4. These results are of particular 
interest to the designer. 

The derivation of the pressure as given in section 1 
is only one of the methods that may be used in this 
problem. As an independent check and as an illustra- 
tion of the application of Busemann’s method of conical 
flow to a nonstationary problem, an alternate derivation 
is given for the pressure distribution in section 4. 

Received August 15, 1948. 
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In a paper by Schwarz! procedures used in oscillating 
airfoil theory are extended to the problem of a wing 
striking a sharp-edged gust at supersonic speed. Re- 
sults for an oscillating down-wash lead to the gust prob- 
lem by a Fourier integral representation. This method 
constitutes a considerable detour and introduces inter- 
mediate results of a transcendental nature which are 
actually not needed and are more complicated than the 
result. It may be verified that the expression derived 
in the present paper for the pressure distribution is 
equivalent to that derived by Schwarz. He does not, 
however, discuss the physical aspects of the problem or 
derive expressions for lift and moment. 


(2 DERIVATION OF THE PRESSURE DISTRIBUTION 


The wing of chord / enters a uniform gust of upward 
velocity v at the supersonic velocity V (Fig. 1). The 
velocity component normal to the wing must remain 
zero, and this condition is equivalent to the generation 
of a velocity normal to the wing which cancels the gust 
velocity (Fig. 2). This may also be considered as a 
“reflection” of the gust on the wing. Because the 
velocity is supersonic, the pressure distribution on one 
side does not influence the pressure on the other, and 
therefore we need only consider the bottom side. The 
pressure distribution on top will be the same except for a 
reversal of sign. For the same reason the pressure dis- 
tribution is not influenced by the trailing edge, and 
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Limits of Integration for Integral (2.6) 





everything is the same as if the wing were of infinite 
chord. The problem is then reduced to finding the 
pressure in a fluid, originally at rest, due to the presence 
of a uniform distribution of normal velocity vp along the 
chord from a point 0, corresponding to the edge of the 
gust, to the leading edge A. If the wing enters the 
gust at time zero the length OA = Vt. 

Now, such a distribution of normal velocity may be 
represented by a distribution of source singularities 
generated continuously at the leading edge and remain- 
ing constant thereafter. The velocity potential of such 
a source appearing suddenly with a constant intensity 
at time ¢, and location x; is 


o= 7 log r — : log E —t)+ V c(t — th)? may A] 


(2.1) 
with 
= (x — x)? + y’ 
This expression satisfies the wave equation 
MER (2.2) 








Ox? = Oy? cr 


Moreover, for sufficiently 
t,), it represents a 


where c = velocity of sound. 
small values of r with respect to c(t — 
velocity field 


o = (vo/m) log r (2.3) 


identical with that of a steady source in the incompres- 
sible flow. Hence, by analogy, it may be concluded that 
the uniform distribution of such sources will produce a 
uniform normal velocity component vp 
The pressure generated by this source along the x axis 
(y = 0) is 
oe pcvg 


po ~es = (2.4 
™ Vet — t)? — (x — x)? ' 





where p = fluid density. Note that the source located 
at x, suddenly appears when the leading edge reaches 
that point, i.e., ata time ¢; = —(x,/V). 
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The local lift 2p at a point x and time ¢, due to the 
uniform distribution of such sources from x, = — Vt to 
x, = 0, is given by 


2p = :. _ f, —— a _ 
. vu Vet + (x1/ 1 — (x — x)? 


(2.5) 
or with the change of variables 
x/ct = €, x/cet = &, c/V = sin p 
I , 
(=> = Mach Number 
sin 4 
. 2 di 
2p = — plvo 2 . : = 3 
r ~via V1 + & om)? — ( — &)* 
(2.6) 


In integrating this expression special attention must be 
given to the limits of integration. The function under 
the integral sign must be taken as vanishing when- 
ever the radical isimaginary. The range of integration 
is therefore limited between the two roots of the equa- 
tion, 


(1 + & sin yw)? — (§ — &)? = 0 


These roots are 


tte 
_ 


2.7) 





£ 1 
(1) — tO = S. 
$l ° ’ sl 
S1 


ny ~ J—sing 





a 


These quantities, plotted as functions of &, are repre- 
sented by two straight lines (Fig. 3) which intersect at a 
point of abscissa 


§ = —(1/sin u) 
and ordinate 
f = —(1/sin ») 


The interval of integration is thus limited to the shaded 
area bordered by the two straight lines [Eq. (2.7) ] and 
the axis (&; = 0). We must therefore distinguish be- 
tween two ranges of values of £: for —(1/sin)< §<—1 


9 £,(2) 
2p = ~ plug ss x 
vis &() 


—— 
V(1+ &sin x)? — (— &)? cose 
Hence, for this interval of {—i.e., between the leading 
edge and the abscissa x = —ct—the lift distribution is 
uniform and independent of x. It may be verified that 
it is identical with the lift on a wing in steady flow with 
an angle of attack »/V. This could have been con- 
cluded immediately, since that portion of the wing can- 
not be influenced by the presence of the gust edge be- 
cause of the finite velocity of propagation of a signal 
originating at this edge. 
The other range of integration —1 < & < 1 corre- 
sponds to points located at distances smaller than ct 
from the gust edge. 
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9 0 dé, 
2p = — plvo = —— “= 
re” Ju V1 + &: sin w)? — (E = 6)? (2.9 
2.9) 
2pcvo ” ( &+ sin yu ) 
a coe ; reer 
@ COS pu 1+ &siny 


In this region the lift distribution depends on the time. 
Substituting § = x/ct, 


2 pcvo x + ctsin p 
2p = - cos7! (= = 2.10) 
7 COS UL ct + x sin pu 
The cos~! branch is between zero and zw. It will be 


observed that the lift distribution depends only on 
x/ct—t.e., there is a similarity law and we may draw a 
one parameter family of lift distribution curves with the 
variable x/ct and the Mach Number as a parameter. 
The appearance of the lift distribution is shown in Fig. 
1. We see that the production of lift is not limited to 
the region of the gust. The lift propagates ahead of the 
edge of the gust with the velocity of sound. The lift 
distribution for various Mach Numbers within the re- 
gion affected by the gust edge is shown qualitatively in 
Fig. 5. 


(3) Lirr AND MOMENT 


It is useful for practical purposes to obtain the values 
of total lift and moment as functions of time for a gust 
acting on a wing of chord /. Because of similarity 
properties, such quantities may be expressed by means 
of nondimensional functions of Vt/l only. In com- 
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puting the time history of lift and moment, we must 
distinguish between three phases. 

Phase 1.—Where the trailing edge is still outside the 
region where lift is produced, (V + c)t < 1. 

Phase 2.—Where the trailing edge is in the region 
influenced by the gust edge, (V + c)t > 11> (V — o)t. 

Phase 3. 
influenced by the gust edge, / < (V — c)t. 


-Where the entire wing is outside the region 


Phase 1 

We integrate the lift distribution 2 over the chord / 
and split the integration into two: intervals, one in 
which the lift distribution is uniform and the other in 


which it is not. Thus the total lift is 


—ct +ct 
baa f pac +2 fi p dx (3.1) 
—Vt —ct 
—1 
2a fi 
—1/sin uw 


p dé + 2ct[pé|ti — 


+1 dp 
2ct E—@e (3.2) 
r f é de 3 


This form is readily integrable. We find 
L= Qpcvol ( Vt/L) (3.3) 


Similarly, the pitching moment about the leading edge 
M, is 


—ct +ct 
M, = 2 ff (Vitxipdr +2 f (Vt+x)pdx 
—Ve —ct 


(3.4) 


or 


=f 1 
waver fo (oases 
ae 2H 
| (G5 +8) a i 
+1 l 2d 
a? ff é. +2) Pit (3.5) 
_; \sin yp dé 


More often we are interested in the stalling moment 
M,,2 about the mid-chord. We find, after integration, 


Vt ea 2 
M,,, = 5 IL — M, = pc! ] (: — ) (3.6) 


/ 


Phase 2 


Similarly, the lift and moment during Phase 2 are 


—ct 1- ve 
L=2 cede +2 fi p dx 


—ct l- Vi 
M, = of (Vits\pdr +2 f (Vt+x)pdx 
—Vi = ¢f 


(3.8) 


(3.7) 





nust 
the 


r10n 
t. 
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dl 
in 
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—1 
L = 2ct id pdé + 2ct[pel’, — 
—1/sin uw 


z dp 
2c &— dé (3.9) 
a fs dé “s 
—] l 
Mr, = eye f ( ‘ + 
~—1/sing \SIN yp 


£) ode + 


2 7 
| (, +8) »| _ 
sin bu 1 
“s l 2d 
cay f i +2) P it (3.10) 
_; \sin p dé 


with s = (1/sin w)[(//Vt) — 1]. We find 


r= =| 3 (: s wad )| 

. sin p l als 
ve Piri l l ' | ™| (3.11) 
ee _— sin p V1 ry oO. 


with L, = 2pcwl/cos uj. The sin—! branch is taken be- 
tween —(7/2) and +(2/2). The stalling moment 
about the mid-chord .1/, /2 is given by 


pcul? al l = Vt sin p 24 
| ey in? * (3.12) 
r\ 1 ysin'n — (G5 -1) roma 


At the end of Phase 2,/ = (V — c)t, and we may verify 
that the above formulas yield L = L, = 2pcul/cos pu 
and M,/. = 0. 


Phase 3 


In this phase the lift and moment remain independent 
of time. We find 1/,,2 = Oand L = L, = 2pcvl/cos p, 
which is the lift in steady-state flight of a wing of angle 
of attack u/V. The lift increases all through Phases 1 
and 2 and reaches its maximum in Phase 3. The time 
history of the mid-chord moment requires special atten- 
tion, as shown in the next section. 


(4) MaximuM VALUE OF THE M1ip-CHORD MOMENT 


With a nondimensional time variable r = Vt/l, the 
mid-chord moment during Phase 1 is [cf. Eq. (3.6) ] 


Mi;2 = pcvol?r(1 “ T) (4 it 


This curve is a parabola with a maximum at 7 = '/2. 


The value of the maximum is 
[Mi /2 mer. ” (1/4) pcvl? (4.2) 


It is interesting to note that during the first phase 
the mid-chord moment is independent of the Mach 
Number. 

During Phase 2 the mid-chord moment is as given by 
Eq. (3.12) 


M'? Ee. l-—r T 
~=7(1 — 7) -] sin : + + 
pcvol” us TSiIn yw 2 


2 | “ 
—— l1—7r\- 

isin? uw — 4.3) 
T \ sail ( T ) \ 


The second phase originates at r = 1/(1 + sin w) and 
terminates at r = 1/(1 — sin uw). It is found that in 
this range the moment 1/,,/2 may go through another 
maximum and that this maximum may be larger than 
its value in the first phase. Value of this maximum and 


“the value of 7 at which it occurs are given in Table 1 for 


various Mach Numbers. 


TABLE 1 
Mach Number [M1/2)maoz 

(1/sin pz) T pcvol? 
1 x x 
ee 2.40 0.331 
ms 1.55 0.281 
1.25 1.10 0.255 
4/m 1.00 1/4 


The maximum of 4/;,2 in the second phase is greater 
than the maximum M,/,. = (1/4)pcw/? in the first 
phase, if (1/ sin uw) < (4/r) = 1.27. The value of the 
absolute maximum of the mid-chord moment is plotted 
against the Mach Number in Fig. 6. This maximum 
is independent of the Mach Number if this Mach 
Number is greater than 1.27. 


(5) ALTERNATE DERIVATION OF THE PRESSURE 


DISTRIBUTION BY THE METHOD OF CONICAL FLOW 

The above results may be derived by an entirely 
different procedure. We may compute the velocity 
field due to the gust reflection on the wing. Because 
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we deal with a supersonic wing velocity, the effect of the 
gust is the same as if the chord were infinite. There- 
fore, the principle of similarity applies, and the velocity 
field is similar at all instants except for a scale factor pro- 
portional with time. In other words, the velocity field 


depends only on the variables 


§ = x/ct, n = y/ct (5.1) 


The velocity field is disturbed by the reflection on the 
wing in a region bounded by the straight lines AFF’ and 
the circle FEF’ centered at the gust edge 0 (Fig. 7): 
In the shaded area AFDF’ the field is uniform and 
corresponds to the steady-state motion of a wing with 
constant angle of attack. The transient field where the 
effect of the gust edge is being felt is inside a circle 
of radius ct centered at 0. The field inside this 
circle may be computed by Busemann’s conical flow 
method. 
By the transformation 


= X sin 8, 
s sin 6 


— = dA cos 8, 
X =scos 68, Y = 
The wave equation [Eq. (2.2)] is transformed to La- 
place’s equation in the X, Y, plane 


(0°¢/OX") + (0°¢/0¥*) = 0 (5.3) 


1+ s*)\ — 
¢(5.2) 


Consider now the components of the velocity field 
= 09¢/dx, v= d¢/dy 


They also satisfy Eq. (5.3). It is convenient to intro- 
duce the complex variable Z = X + 7¥Y. Let us in- 
vestigate the velocity field on the bottom side of the 
wing. Thev component of the velocity field is 


y= Re™ [log (Z - Z:) + log (Z — Zz) — log Z] 
(5.4) 


where Re = real part of; Z, = ie“; Z, = —ie~™; sin 
u = c/V. This expression satisfies Eq. (5.3) and the 
boundary value of v on the circle and the wing—i.e., 
—v on ODF’ andv = C on F’EO. In terms of real 
quantities, 


=% J [s — 
- tan~ “Us + (1/ 


9= 


(1/s)] ] sin 6 i] 
s)] cos 0 +2 2 sin P| 


(5.5) 





Now we are interested in the pressure distribution on 
the wing. This pressure distribution for 7 = 0 may be 
derived from the above expression for v by making use 
of the equations of motion and continuity 


p (0u/dt) = —(Op/dx) 


(> ‘. *) i. (5) 
P\dx dv) sc? “\Ot 


By elimination of « and transformation of variables, 
we find that, on the x axis, cos@ = +1. This implies 


(5.6) 


cos 0(0v/00) = [(A? — pc |(Op/Ord) (5.7) 
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By integration 
(5.9) 





pcvo , (! ert sin uM 
= ———= Con 
© COS 1+ ésin pu 


which coincides with expression (2.9) above. 


CONCLUSIONS 


It has been shown that the pressure distribution on a 
supersonic wing striking a sharp-edged gust may be ex- 
pressed by a simple formula. This pressure distribu- 
tion is obtained by direct integration of a variable dis- 
tribution of sources in a fluid at rest. It is also shown 
that the same result is obtainable from Busemann’s 
method of conical flow. From the time history of total 
lift and moment it is concluded that the largest value of 
the total lift is reached in the last phase—i.e., when a 
steady flow has been established—while for the mid- 
chord moment a maximum value [J/1/2],,,.. = (1/4) pcvl* 
is reached if the Mach Number is larger than 4/7. 
This maximum moment is independent of the Mach 
Number. However, for Mach Numbers between 4/7 
and 1, the maximum mid-chord moment varies with the 
Mach Number and becomes infinite at Mach Number 1. 

(Continued on page 310) 
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Column Curves for Magnesium-Alloy Sheet 


EVAN H. SCHUETTE* 


The Dow Chemical Company 


ABSTRACT 
Results are presented for flat-end tests of 129 columns cut from 
flat FS-1h, FS-la, and Ma magnesium-alloy sheet 0.250 in. thick 
and for clamped-end tests of 26 columns cut from flat FS-lh 
sheet 0.064 in. thick. The best empirical fits to the data were 
obtained with the parabolic formula 
P (KF .,)?(L’/R)? 


— = KF, — 
A 4r2E 


with a horizontal cutoff at P/A = Fry. 


INTRODUCTION 


pee BEHAVIOR OF SIMPLE COLUMNS within the elastic 
range is well understood. For predicting column 
strength beyond the elastic range, a number of methods, 
both empirical and theoretical, have been proposed. 
In any case, it is desirable to have test data to substan- 
tiate whatever method is selected. Reference 1 pre- 
sented such test data and empirical formulas for mag- 
nesium-alloy extrusions. The present report gives test 
data for 155 columns of magnesium-alloy sheet and 
shows the correlation with the tangent modulus column 
curve and empirical parabolic curves. 


SYMBOLS 


= area of column cross section, sq.in. 
E = modulus of elasticity, ksi 

E, tangent modulus, ksi 

F., = compressive yield stress, ksi 

K 


= empirical constant 

= column length, in. 

= column effective pin-ended length, in. 

= load, kips 

R = minimum radius of gyration of column cross section, in. 

fixity coefficient in the Euler column formula when given 
in the form P/A = cr?E/(L/R)? 

n = empirical constant 

ksi = kips (thousands of pounds) per square inch 


SNS 


ll 


SPECIMENS AND METHOD OF TESTING 


The materials tested were FS-lh, FS-la, and Ma. 
Chemical compositions are as follows: 


‘ 


Other 
% % vi % “% Impur- 
// Mn, a Si, Cu Ni, Fe, ities, 
Alloy Al Min Zn Max. Max. Max. Max. Max. Mg 
FS-1 2.5- 0.20 0.7- 0.3 0.05 0.005 0.005 0.3 Remainder 
3.5 
0.3 0.05 0.01 ina 0.3 Remainder 


M - .20 


The suffix ‘“h’’ denotes the hard condition; ‘‘a,’’ the 


annealed condition. Applicable specifications and 


Received December 20, 1948. 
* Division of Tests, Magnesium Laboratories. 


specification properties for the materials are listed 


below: 
Min Min. Com 
Min. Tensile pression 
Tensile Yield Min Vield 
Specifica- Strength Strength Per Cent Strength* 
Material tion (ksi) (ksi) Elongation (ksi) 
FS-1h AN-M-29 39 29 4.0 24 
FS-la AN-M-29 32 “— 12.0 ‘ 
Ma AN-M-30 28 s 12.0 





* Not specified. 


Two sets of test data were taken. A brief series of 26 
clamped-end tests of FS-1h were run in accordance with 
the method outlined in reference 2. The specimens for 
these tests were nominally 0.064 in. thick and 0.400 in. 
wide. The tests were made on a 60,000-lb.-capacity 
hydraulic compression machine using the 2,400-lb. 
scale, which could be read to the nearest pound. This 
is equivalent to an increment of about 0.04 ksi on the 
column. Specimen lengths and test results are given 
in Table 1. 

The remaining 129 tests, covering FS-lh, FS-la, and 
Ma sheet, were on specimens nominally 0.25 in. thick 
and 1.00 in. wide. Ends and edges were milled flat, 
square, and parallel (see reference 1 for comments on 
the machining accuracy). The tests were made on the 
same testing machine as the previous series, with the 
difference that all columns taking a load of more than 
2,000 Ibs. were tested on the 12,000-lb. range. The 
columns were tested with their upper ends bearing 
directly against the upper head of the testing machine. 
Because the upper head was not adjustable, and ram 
stroke was limited, length adjustment was made by 
placing a series of machine-fitted blocks on top of the 
ram, with the lower end of the specimen bearing against 
the uppermost of these. Some added degree of testing 
machine flexibility, or ‘‘slop,’’ was doubtless introduced 
by use of these blocks; its effect, however, could not 
readily be calculated. Specimen lengths and test re- 
sults for this series of tests are given in Table 2. 

An effort was made to obtain tangent-modulus curves 
applicable to the rate of loading used for column test- 
ing, by taking them from autographically recorded 
stress-strain curves. The curves were not steady 
enough, however, to permit making accurate slope 
determinations. Consequently, while these curves suf- 
ficed to prescribe compressive yield stresses, the tan- 
gent-modulus values used had to be obtained from 
typical stress-strain curves for the materials.* The 
latter were not autographically recorded and, conse- 
quently, correspond to a slower rate of loading than 
that used for the column tests. 
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TABLE 1 
Clamped-End Test Results 
Fy Length P/A P/A V Fey(L’/R 

Alloy Direction (ksi) (in.) (ksi) L’'/R Fey (ksi)®.5 

FS-1lh Long. 26.6 0.60 26.60 15.8 1.000 81.5 
0.75 25.55 19.9 0.960 102.6 
0.95 25.05 25.4 0.941 131.0 
1.10 25.00 31.1 0.940 160.4 
1.30 23.20 35.2 0.872 181.8 
1.50 22.05 42.8 0.830 220.6 
1.70 20.42 46.3 0.768 239.0 
1.85 19.60 §2.5 0.737 270.8 
2.20 15.23 59.5 0.573 307.0 
2.58 12.30 70.4 0.462 363.0 
3.70 6.56 98.4 0.247 507 .0 
4.50 4.08 126.0 0.1538 650.0 

FS-1h Trans. 28.3 0.60 29.50 15.9 1.042 84.6 
0.75 27.40 21.2 0.969 112.7 
0.95 26.30 25.6 0.930 136.2 
1.10 26.00 29.9 0.919 159.0 
1.30 26.00 37.6 0.919 200.0 
1.50 22.60 43.9 0.799 233.6 
1.70 21.20 50.1 0.750 266.6 
1.85 18.40 53.9 0.650 286 .6 
2.20 14.22 60.3 0.503 321.0 
2.58 12.40 70.0 0.438 372.0 
3.00 9.70 81.0 0.343 431.0 
3.35 8.00 89.0 0.283 473.0 
3.70 5.56 107.5 0.197 572.0 
4.45 4.03 125.0 0.143 665.0 


METHOD OF REDUCING DaTA 


Results of the 26 clamped-end tests were reduced to 
pin-ended values according to the procedure of refer- 
ence 2. For the flat-end tests, an average value of 
(L’/R)/(L/t) was obtained from the test results that 
fell within the elastic range. This value was then used 
in reducing all the data. The corresponding average 
values of end fixity coefficient c were 3.72 for the 
clamped-end tests and 3.46 for the flat-end tests. These 
values appear reasonable, considering the out-of-flat- 
ness that is likely in columns cut from nominally flat 


sheet. 


CURVE FITTING 


Considerable success was obtained in reference 1 in 
using generalized parameters for plotting column data. 
The failing stress was given in terms of the compression 
yield stress; that is, as (P/A)/F,,. By dividing both 
sides of the Euler formula by F,,, 


(P/A)/Fey = wE/Fe(L'/R)? 


the second parameter, F,,(L’/R)? or v ‘Fey(L’ R), is 
obtained. These same parameters were used in the 
present investigation. 
for fit, the parameter V F.,(L’ R) was raised to suitable 
powers, so that the curve to be checked could be drawn 
as a straight line; in some cases the inverse form, also, 
F.,/(P/A), was used to simplify fitting. The following 
types of formula were tried. 


In checking the several formulas 


Hyperbolic: 
(P/A)/Fy = K/(V F.,(L’/R)\" 
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Rankine: 


(P/A)/Fey = Ki/{1 + Ko[W F.,(L’/R)"} 
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Alloy 


FS-1lh 


FS-1h 


FS-la 


COLUMN CURVES 


Direction 


Long. 


Trans, 


Long. 


28.4 


16,5 





TABLE 2 
Flat-End Test Results 


Length P/A 
(in.) (ksi) 
2.00 25.80 
2.00 25.80 
2.77 25.52 
2.77 25.70 
3.72 23.64 
3.72 25.54 
4.85 23.20 
4.85 25.38 
5.45 25.16 
5.44 23.36 
6.13 23.64 
6.13 23.80 
6.76 21.30 
6.76 21.30 
7.52 18.08 
7.52 19.95 
8.18 16.83 
8.19 16.49 
8.95 14.12 
8.95 15.30 
9.64 13.39 
9.64 13.70 
10.37 11.09 
10.37 11.54 
11.69 8.45 
11.68 9.56 
13.01 7.53 
13.01 7.61 
14.52 5.59 
14.52 5.88 
16.19 4.56 
16,19 4.37 
2.00 30.20 
2.78 27.46 
3.72 27.00 
4.85 25.36 
5.45 26.40 
9.13 22.72 
5.76 22.02 
7.53 18.03 
8.19 16.98 
8.95 14.20 
9.64 14.00 
10.37 11.61 
11.68 9.37 
13.01 7.87 
14.52 5.57 
16.19 4.32 
2.74 16.60 
2.76 16.92 
3.81 16.57 
3.82 iy a7 
5.08 16.86 
5.08 16.66 
6.55 15.90 
6.54 iy. 
0.731 13.73 
0.731 14.77 
8.19 13.40 
8.19 16.61 
0.904 13.42 
0.904 12.78 
9.96 11.18 
9.96 12.79 
10.82 10.22 
10.82 10.18 
11.74 9.84 
11.74 9.38 
12.60 7.61 
12.60 7.75 
13.49 6.86 
13.49 7.05 
15.04 5.72 
15.04 5.13 
16.81 3.92 
16.81 4.57 


(Table 2 continued on page 304) 
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TABLE 2 (Continued) 


Fey Length 
Alloy Direction (ksi) (in.) 
FS-la Trans. 17.4 2.76 
3.81 
5.08 
6.55 
7.31 
8.19 
9.04 
9.96 
10.81 
11.74 
12.59 
13.49 
15.04 
16.81 
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Parabolic (including straight line): 

{(P/A)/Fe = Ki — Ko[V F,(L'/R)]"} 
Tangent modulus: 

(P/A)/Fy = 2°E,/F.(L'/R)? 


Of these, the parabolic (with » = 2, A, = 1.00 to 
1.20, Ay adjusted to give tangency with the Euler 
curve) and the tangent modulus formulas gave the best 
fits. These are shown in Figs. 1-3. In Fig. 3, it will be 
noted that the tangent-modulus curve is unduly con- 
servative for Ma. This failure to get good agreement, 
however, may have been due to differences between the 
material used for column tests and that used for deter- 
mining the tangent modulus. Moreover, it is known 
that, of the three materials, Ma shows the greatest 


V Fe,(L’/R 


P/A P/A 

(ksi) L/h Fey (ksi)°-5 
17.85 20.1 1.026 83.9 
17.70 27.8 1.018 116.0 
17.73 36.9 1.019 153.9 
17.98 47.9 1.033 199.7 
16.12 53.4 0.926 222.6 
17.08 59.5 0.981 248 .2 
13.36 66.0 0.767 275.4 
12.10 72.6 0.695 302.6 
9.94 79.0 0.571 329.4 
8.22 85.6 0,472 357.0 
8.41 91.7 0.484 382.2 
6.86 98.2 0.394 409.0 
5.32 109.6 0.306 457.0 
3.99 123.1 0.229 514.0 
11.96 26.2 0.980 91.4 
12.22 26.2 1.002 91.4 
11.32 35.6 0.928 124.2 
11.73 35.6 0.961 124.2 
11.52 46.8 0.945 163.2 
12.15 46.8 0.995 163.2 
10.76 53.0 0, 882 185.0 
10.38 52.9 0.851 184.7 
11.63 59.9 0.953 209.0 
12.02 59.6 0.985 208.0 
9.39 66.5 0.770 232.0 
9.10 66.5 0.746 232.0 
8.96 74.1 0.734 258.6 
9.66 74.1 0.791 258.6 
9.06 81.0 0.743 282.6 
8.87 81.0 0.727 282.6 
7.86 89.0 0.644 310.6 
7.70 89.0 0.631 310.6 
5.57 98.5 0.457 343.8 
6.94 96.3 0.569 336.0 
5.69 104.9 0.466 366.0 
5.72 104.4 0.469 364.2 
4.11 118.8 0.337 415.0 
4.21 119.3 0.345 417.0 
3.46 134.2 0.284 468 .0 
3.00 133.8 0,246 467.0 
11.34 26.2 0.970 89.5 
11.27 35.5 0.963 121.4 
10.46 46.9 0.894 160.3 
11.38 53.2 0.973 182.0 
10.04 60.9 0.859 208.0 
9.96 67.2 0.851 230.0 
8.24 75.1 0.704 256.8 
8.62 81.6 0.737 279.2 
6.36 90.2 0.543 308.4 
6.65 96.3 0.569 329.4 
5.01 104.4 0.428 357.2 
4.85 119.2 0.414 408.0 
3.36 132.7 0.287 454.0 


decrease in tensile yield strength with decrease in rate 
of loading. It is possible, therefore, that compressive 
properties are similarly affected and that a tangent- 
modulus column curve based on a rate of loading cor- 
responding to that used in the column tests would pro- 
vide an improved fit to the data. 

In all cases, a large amount of scatter is evident, so 
that neither type of curve gives an accurate predic- 
tion for all the test results. Conversely, then, it is 
unimportant which of the two curves is selected for 
use. 


MINIMUM PROPERTY CURVES 


The generalized form in which the column formulas 
are written makes possible a direct adjustment for any 
specified values of compressive yield stress. Fig. 4 
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FIG. 2- COLUMN-TEST RESULTS FOR FS-la SHEET 


shows such adjusted curves, in the conventional P/A 
vs. L’/R plot, for FS-1h and FS-1Wh. The latter is a 
calcium-free, weldable modification of FS-1h. 


CONCLUSIONS 


The best purely empirical short-column formula indi- 
cated by the test data is the parabolic formula 


. kh vant ho R)? 
A Io? 
with a horizontal cutoff at P/A = F,,. Values of K for 
the several alloys are as follows: 
FS-1ih longitudinal K = 1.05 
FS-1h transverse 1.05 
FS-la longitudinal 1.10 
FS-la transverse 1.20 
Ma longitudinal 1.00 
Ma _ transverse 1.00 


The tangent-modulus column curve also provides 
good fits to test data for FS-1h and FS-la but is overly 
conservative when applied to Ma. The failure to get 
good agreement for Ma may be due to differences be- 
tween the material used for column tests and that used 
in determining the tangent modulus or to differences 
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FIG. 3— COLUMN-TEST RESULTS FOR Ma SHEET 
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FIG. 4- MINIMUM-PROPERTY COLUMN CURVES FOR FS-th @ FS-lWh 


in rate of loading between the column tests and the 
stress-strain curve determinations. 

All the test results were characterized by large scatter, 
probably a result of random out-of-flatness of the 
columns. 


(Continued on page 310) 








The Potential Flow Around Given Wing 


Sections 


KWEI-LIEN CHANG* 
Bureau of Aircraft Industry, China 


ABSTRACT 


By the use of Wilton-Morris’ parametric representation! and 
the K&arman-Trefftz transformation,? the complex potential 
around a given wing section can be given as an exponential series. 
This series is exact in the sense that the actual profile instead of 
the approximate one is used for its developing. 


INTRODUCTION 


= SO-CALLED DIRECT PROBLEM—i.e., to find the 
potential flow around a given wing section 

been solved by various methods, notably Karman- 
Trefftz approximate method? and Theodorsen’s e-func- 
tion method.* The Karman-Trefftz method involves 
an approximation in using g instead of log, (1 + g) and 
it becomes heavily involved when terms with higher 
Theodorsen’s 


has 


powers are retained as by Hohndorf.* 
method gives solution on the boundary only; its calcu- 
lation depends on graphical or numerical differentiation 
that is rather doubtful in accuracy. In addition, the 
accuracy near the trailing edge is further affected by the 
preliminary use of the Joukowski transformation instead 
of the Karman-Trefftz transformation. The method 
presented below gives accurate results near the trailing 
edge, as well as near the leading edge, and it applies 
to the field as well as on the boundary. The labor in- 
volved is just equal to that required in the first ap- 
proximation of the Karman-Trefftz method, and the 
accuracy achieved depends on the number of terms 
retained in the series, of which usually a few terms give 


accurate enough results. 


NOTATIONS (REFER TO FIG. 1) 


x,y = rectangular coordinates in the s-plane of the physical 
profile 

r,@ = polar coordinates in the z-plane 

X, Y = rectangular coordinates in the Z-plane of the circular- 
like profile after the Karman-Trefftz transforma- 
tion 

p,@ = polar coordinates in the Z-plane 

ft, = rectangular coordinates in the ¢-plane 

Ww = complex potential 

r = circulation 

V = velocity at infinity 

a = the angle of attack 

c = chord 

T = trailing edge angle in radians 

R = radius of curvature at the leading edge 

k = [2 — (r/r)],a constant in Karman-Trefftz transforma- 
tion 


Received September 8, 1948. 
* Aerodynamics Engineer. 


306 


l = a linear dimension in Karman-Trefftz transformation 
so selected as to adjust the trailing edge of the pro- 
file at (—&l, 0) 
Be = Conjugate function of the function B, 
de, @1, A2..., dy, be..., An, Bn, A,’, Bn’ equal real constants, 
Subscript ‘‘zero”’ signifies points on the boundary. 


n 


(1) INriTIAL TRANSFORMATION 


This initial step is required to render the method of 
practical use; otherwise the number of terms required 
in the series may be so many as to make the method 
useless. 

Let the points in the physical plane be denoted by 
z= x-+iy. The transformation 


(gs — kl)/(2+ kl) = [((Z-D)/(Z+D)]* (i) 


transforms the given wing section in the z-plane into 
a circular-like profile in the Z-plane.” & is related to the 
trailing edge angle 7 by 


k = 2 — (r/m) 


and / and the origin are so selected as to make (—8i, 
0) coincide with the trailing edge, and (+i, 0) is then 
near to the point of maximum curvature at the leading 
edge but is just inside the profile. 
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POTENTIAL FLOW 


Because of accuracy required, the geometrical method 
is not precise enough, so that analytical formulas in 
rectangular coordinates are derived in the following 


AROUND 
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In Eqs. (4) and (5) the negative square root is discarded 
because it transforms the region at infinity in the z-plane 
to the origin in the Z-plane. 


expressions, since the conventional airfoil is usually 


expressed by its ordinates along the chord. 


Eq. (1) is equivalent to the following two equations 
giving the magnitude and 


[Eqs. (2) and (3)], Eq. (2) 
Eq. (3) giving the directions. 
e~oar+7* _ E —1)?+ =| 
(x + kl)? + y? ~ L(¥ + 22+ ¥? 


tan= ( y =) — tan i( z ) = 
x- 


\x + i 


Y Y 
&| tan : (. £22 ) — tan (- + | (3) 


By Eqs. (2) and (3) and the substitutions 


(2) PARAMETRIC REPRESENTATIONS BY FOURIER 
SERIES 


The polar coordinates in the Z-plane are given by 
(2) p = (X?+ Y*)’ 
@ = tan= (¥/X) 


By the points (x0, yo) on the given section, po can be 
plotted against ¢» from 0 to 27 and expressed as a 
Fourier series, 


po = A + pa a, COS Ndy + > b, sin n¢do (6) 
I I 


2/k cas )2 2 | /k 
k,?2 = (“) on E kl)’ + ¥ | where a's and 0’s are real constants. Since, 
‘< (x + Rl)? + y? 
X = pcos ¢ 
: ie i Si x — ki Y = psin @ 
Vv 
tz —1 : So, 
1 (. + a) | 
ay Qe = “J ane +1 
C=1- (1 —k2)/0 +k)? me oe («5 + *) cos bo + x 
Co = tan? ky + [(1 — k,?)/(1 + R,?)]? be an + be 
we obtain, cos do + > sin do + :- : ——— sin ndo 
, 1— k\01 + V1 + (e/c,) b, be has ~ &- 
X = hs dines . Rewe=-4~ 60 Sit, BeeB st 
4 i(; + “A | o/s Tay Bests ere 2 . 
tan ke (1 + V1 + (o/c, P “1 Gs 
fat n+ + (ale 2 (5) cos Ngo + («. -< *) in go + he : toner sin 2do 
Ce/ Cr “ 
With the same notation as in reference 1, put 
é = — go 
We obtain, 
Lo = Xo + 1Vo 
1 lo - n + aAn-— by ‘ : b, T b . 
= E + (a. + *) cosé+ >> noon a cos n~i — ; siné — >> at = . * sin nt | aa 


b by bn oi 2 as — An,-1 . 
i| ** cos —§ + wf, Sten cos né — (« _ *) sin e+ > pent. ~ at sin ne | (7) 


9 9 


By the two identities below, 


(A, + iB,)e~™ = (A, cos nt + B, sin nt) + i(B,cos nt — A, sin né) 


(A,,’ 4 iB, Netting — (A,’ 


cos n& — B,’ sin nt) + i(B,' cos né + A,’ 


sin n£) 


we can write Eq. (7) as Eq. (8) in the exponential form 


LZ = fie®) 


b . bo o |] a l i 
= (2: + i?) + age~* + (S + i) wr 4 2 («. 1 — ib,- :) or + 2s («. ti + ib. +1) e**= (8) 
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(3) THe COMPLEX POTENTIAL 

After Z, is expressed as a function of & exponentially 
to represent the profile, the complex potential can be 
written out by the idea of parametric mapping.’ Sub- 
stituting ¢ = & + im for &in Eq. (8), we obtain a func- 
tion 


, “it ay . by , 
Z= fle") =\|— +3 + ae 
bs " | » 

7 Je ds An-1 — 1b,-1) e~*™ + 


| («. 11 + i, -) ett (9) 


5 » 


which maps the region outside the circular-like profile 
in the infinite strip0 < & < 27,0 <» < © inthe ¢-plane. 

The required complex potential in the Z-plane must 
satisfy the following two conditions in addition to its 
analytical property: 


(a) The boundary of the circular-like profile—i.e., 
» = 0—tust be a stream line. 
(b) dW/dZ = (dW/d¢)(dt/dZ) must be equal to the 


negative conjugate of the free-stream velocity, + Ve™, 
when Z = ~©—i.e., when 7 = o, 


Following reference 1, write 
f(e®) = By(e*®) + Ba(e*) 


where B,(e*) represents that part of f(e") with positive 
exponents and B.(e*) that part with negative exponents. 
Then the complex potential can be shown as: 


W = + Ver*Bi(e®) + Ve-B(e-*) + Te (10) 


by following the same line of derivation as in reference 
1 for moving axes without circulation. After substitut- 
ing B,(e*) and B.(e-*) from Eq. (9), Eq. (10) becomes, 
W = + Vet™ [ace is + > (1 2(Gy- = 


ib, —1)e—™] + Ve-[acet® + > (1/2)(an-1 + 


ib,—:yet™] + TE (11) 


in which the constant term has been neglected because it 
does not influence the flow pattern. 

Finally, we come to the circulation term Tf. It is 
shown in reference 1 that 'f = T'é + 2In gives the term 
for circulation around the profile, because the imaginary 
part vanishes on the profile and the real part decreases 
by 27I when we go once around the profile in the posi- 
tive direction. The magnitude of I will be determined 
by the Kutta condition at the trailing edge—i.e., 


dW/dZ = (dW /dgé)(dé/dZ) = 0 


atyn = 0, = 7; or 
dW/d§ = 0 


S 
See 
Il 


7; since 


dg/dZ # 0 


atyn = 
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atn = 0, = zw. So we obtain 


= & 
wie Se 


‘e “4 we ive | ave + 
n 1 4 n 

Us («. i+ ibaa) ev | (12) 
— 9 


PROCEDURE OF CALCULATION AND DISCUSSIONS 


l= +iVer"| ae . 


(0. -1 —1b, 


(4) 


The procedure for finding the complex potential of 
a given witig section can be summarized as: 

(a) By Egs. (4) and (5) to transform the given sec- 
tion into a circular-like profile. 

(b) After changing the coordinates Xo and Y, to 
po and go, po is expressed in a Fourier series of ¢ as 
Eq. (6). 

(c) With coefficients of Eq. (6), the complex poten- 
tial of the circular-like profile is directly written out by 
Eqs. (11) and (12), and then the complex potential of 
the given section can be obtained by substitutions of 
¢, Z according to Eqs. (9) and (1), which are usually 
not necessary to carry out. 

The most time-consuming work is involved in step 
(a); step (b) can be carried out quickly and accurately 
by using the established scheme for 24, 36, or 48 ordi- 
nates; and the step (c) is just some simple substitutions. 
The accuracy achieved depends on how many terms are 
used in step (b); as illustrated in next paragraph, for 
that symmetrical section, five or six terms in the Fourier 
series will give sufficiently accurate results, about 0.05 
per cent. 

It can be seen that the labor involved is just equal 
to that required in the Karman-Trefftz method for the 
first approximation; when the region near the trailing 
edge is not considered important, the Joukowski trans- 
formation may be used as an initial step to reduce the 
labor. 


(5) NUMERICAL EXAMPLE 


Step (a)—The given section, as presented in columns 
(1) and (2), Table 1, has its trailing edge angle, 
7 = 10.8° = 0.06 rad. 


1.94 


k= 2 —(r/r) = 


Now we intend to have the trailing edge represented by 
(—kil, 0) and have (+i, 0) to represent a point near 
the leading edge. Theoretically, for k = 2, it is best to 
have (ki, 0) to represent the point halfway between the 
leading edge and its center of curvature ;’ for 1 < k < 2, 
(+l, 0) should be chosen to lie between the halfway 
point and the leading edge to give the best results. 
For the present case, the value of & is so near to 2 that 
the halfway point will be used. Since, R = 0.686 per 
cent chord, 
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2kl = [100 — (1/2) 0.686]% c = 99.657% c 
1 = 25.685% c 
With this selection of k/, the coordinates of the given 
section are presented in columns (3) and (4) referred 
to the new origin. By Eqs. (4) and (5) the transformed 
circular-like profile is calculated, and it is presented in 
columns (5) and (6) in polar coordinates. 

Step (b)—Plotting the values in column (5) against 
those in column (6), as shown in Fig. 2, then, by using a 
48-ordinate scheme, pp is expressed as a Fourier series 
of @» as follows: 


i ii 
j = ao + > an cos nd + D> dy sin nd 
1 i 


where 

a, = +0.05863 
a3 = +0.00291 
as = +0.00093 


a = +1.05183 
a2 +0.00935 
a, = +0.00168 
+0.00045 


ae 


and by, b:, ...b, = 0. After checking some points be- 
tween po//, represented by Fourier series on to sixth 
harmonics, with that in columns (5) and (6), it is found 
that the error involved in the Fourier representation is 
only an order of 0.05 per cent. 

Step (c)—Finally, the complex potential is written 
out by Eqs. (11) and (12) as: 


W = + Vet[1.05183e—* + 0.02932e-* +... J+ 
Ve-*[1.05183e+® + 0.02932e+# +...J0 + Te 





TABLE | 
Given Symmetric Transformed 
Section with Given Section Profile 
Origin at with Selected in Polar 
—Leading Edge— —_——Origin——-—~ —Coordinates— 
Distance 
from e 
leading 
edge in Ordinate Xp in Yo in ¢ in 
%C mAc » ET %C P./l degrees 
0 0 —50.171 1. 1260 0 
0.5 0.8254 —49.671 1.1245 8.10 
0.75 1.0092 —49.421 1.1234 9.75 
1.25 1.2985 —48.921 1.1215 12.90 
2.50 1.8205 —47.671 1.1170 18.45 
5.0 2.5293 —45.171 1.1090 26.55 
7.5 3.0412 — 42.671 1.1010 32.80 
10.0 3.4450 —40.171 1.0975 36.85 
15 4.0505 —35.171 1.0885 45.50 
20 4.4727 —30.171 1.0795 53.60 
25 4.7595 —25.171 1.0730 60.15 
30 4.9329 —20.171 Same 1.0665 66.30 
35 4.9992 —15.171 as 1.0590 73.65 
40 4.9544 —10.171 Column 1.0535 79.93 
45 4.7699 — 5.171 2 1.0470 85.83 
50 4.4920 — 0.171 1.0420 91.92 
55 4.1465 — 4.829 1.0390 95.20 
60 3.7507 — 9.829 1.0335 101.63 
65 3.3179 — 14.829 1.0285 107.20 
70 2.8598 —19.829 1.0240 112.70 
75 2.3873 — 24.829 1.0195 119.60 
80 1.9102 — 29.829 1.0150 127.50 
85 1.4326 —34.829 1.0115 133.33 
90 0.9551 —39.829 1.0075 142.95 
95 0.4775 —44.829 1.0040 153.10 
100 0 —49_ 829 1.0000 180,00 
rt = 10.8°, ki = 49.829% C | = 25.685% C 


R = 0.686% C 


Note: The numerical figures in columns (5) and (6) are not 
considered accurate because the Karman-Trefitz transformation 
was carried out with improper numerical tables; nevertheless, it 
can still be used to demonstrate the main aspects of the method. 
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ABSTRACT 


The usual approximations made in calculating the pressure 
changes of a compressible fluid flowing through heat-exchanger 
passages introduce appreciable errors in the range of high flow 
Mach Numbers and high rates of heating. Existing methods for 
obtaining accurate results over the entire Mach Number and 
rate of heating range require numerical integration for each set 
of conditions and are thus too tedious for general application. 

In the present paper an analysis is made of the compressible- 
flow variations occurring in heat-exchanger passages. The re- 
sults of the analysis describe the flow and heating characteristics 
for which specific flow passages can be treated as segments of 
generalized flow systems. The graphical representation of the 
flow variations in the generalized flow systems can then be uti 
lized as working charts to determine directly the pressure changes 
occurring in any given flow passage for any given flow and heating 
conditions. Working charts that have been constructed on the 
basis of these results to handle the case of air heated at constant 
wall temperature under turbulent fluid-flow conditions are dis- 
cussed. A method is given of incorporating the effect on the 
heat-exchanger flow process of high temperature differential 
between passage-wall and fluid as based on recent N.A.C.A. 
experimental data. Good agreement is obtained between the 
experimental and the chart pressure-drop values for passage- 
wall average temperatures as high as 1,700°R. (experimental 
limit) and for flow Mach Numbers ranging from 0.32 to 1.00 
choke) at the passage exit. 


INTRODUCTION 


—~ RATIONAL DESIGN and calculation of perform- 
ance of aircraft heat exchangers, wherein heat is 
added to a high-speed compressible fluid stream, re- 
quires, first, the knowledge of the appropriate heat- 
transfer and friction-drag coefficients and second, the 
mathematical means of utilizing these coefficients to 
describe accurately the pressure and temperature varia- 
tions of the fluid along the length of the heat-exchanger 
flow passage. 

One form of the differential momentum equation de- 
scribing the one-dimensional steady-state motion of a 
compressible fluid in a constant-area flow passage under 
the simultaneous influence of friction and heat addition 
is 


d(mV + pA) + dD; = 0 (1) 


Presented at the Heat Transfer and Fluid Mechanics Institute, 
Los Angeles, June 21-23, 1948. Received September 1, 1948. 
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where m is the mass flow; V, the velocity; /, the static 
pressure; A, the cross-sectional flow area; and D,, 
the drag force due to friction. The grouping (mV + 
pA) is referred to as the total momentum of the flow 
(m\ is the momentum force, and pA is the pressure 
force). In Eq. (1) the velocity I’, static-pressure p, 
and drag force D, are functions of distance along the 
flow passage. 

No exact closed-form integration of this differential 
equation has yet been obtained, even for specially 
chosen boundary conditions. Two general methods of 
approach have been used for obtaining solution of the 
differential equation. 

The first method has been to introduce approxima- 
tions (as, for example, assuming constant average fluid 
density in the flow passage) which permit formal inte- 
gration of Eq. (1) and yet which introduce little error 
in some restricted range of interest. This method has 
been successfully applied in the range of low flow Mach 
Numbers and low rates of heat input to the fluid for 
which conditions the flow variations along the heat- 
exchanger passage are not critically influenced by com- 
pressibility effects. No satisfactory application of this 
method of solution has been possible for the high-speed 
high-heating range. 

The second method has been to integrate numerically 
Eq. (1) for each specific condition of heat-exchanger 
operation. Accurate results are thus obtained over the 
entire flow and heating range. This method of solution 
is, however, much too tedious, particularly when in- 
vestigation of a large number of heat-exchanger oper- 
ating conditions is required. 

The problem at hand is, therefore, to provide a con- 
venient means for accurately determining the pressure 
variations of a compressible fluid flowing in a heat- 
exchanger passage under conditions of high flow Mach 
Numbers and high rates of heating. 

This problem has been analyzed and the results ap- 
plied to the specific case of air flowing in turbulent 
motion through smooth-walled passages and heated at 
constant wall temperature conditions (Valerino'). A 
general treatment of the analysis and a discussion of the 
physical significance of the results are presented herein. 
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The analysis results are compared with experimental 


data. 
DISCUSSION 
Eq. (2) is a form of the conservation of energy equa- 
tion 
t(/T =1— [(y — 1)/2y](V/V gRT)? (2) 
where 7 is the absolute total-temperature; ¢, the abso- 
lute static-temperature; y, the ratio of the specific 
heats; g, the mass conversion factor; and R, the gas 
constant. 
Eq. (3) is the isentropic relation 
P/p = (T/)Y (3 
where P is the total pressure. 
Eq. (4) is a statement of the conservation of mass 
principle when introducing the perfect gas law: 
t/T = V gRT) 


Eq. (5) can then be written from Eqs. (2) and (4): 


(pA ‘mvV/gRT)(V (4) 


mV + pA 14 [(y +1)/2)(V/-WgRT)?] (5 
eae — = — ») 
mv/gRT V/V gRT 


The parameter (mV + pA)/m+/gRT is referred to 
as the total-momentum parameter of the flow. Eqs. 
(2) through (5) show that, when the total-momentum 
parameter (mV + pA), mV/gRT and the total-tem- 
perature 7 of a fluid having a certain y are known, all 
the fluid-flow conditions (p, P, V, and ¢) are uniquely 
specified. The problem of describing the flow process 
in a heat-exchanger passage thus reduces to one of de- 


(” V + E) Y 
= on 
. mv/gRT _ 4 


dT 7 


If Reynolds analogy is assumed to hold, (h/c,pgV) + 
(F/2) is equal to a constant; hence, Eq. (9) involves 
only the variables (mV + pA)/mvV/gRT, T, Tw, and y. 

The conclusion can then be drawn that, for a fluid 
having a certain value of y, the variation of (mV + 
pA)/m+/gRT during the flow process is uniquely deter- 
mined by the variation of T during the process, pro- 
vided that 7,, is constant or is expressible as a function 
of only (mV + pA)/mV/gRT and T. 

Thus, with appropriate restrictions on the passage- 
wall temperature, the variations of the fluid-flow condi- 
tions between any two stations in a flow system are 
uniquely determined by the flow conditions at one of 
the stations and by the variation of total temperature 
between the two stations and, in particular, are explic- 
itly independent of the absolute positions of the two 
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fining the variations of the total-momentum parameter 

and the total temperature along the flow passage. 
The variation of the total-momentum parameter 

(mV + pA) ‘mV gRT is due to the variation of both 

the total momentum mV + pA and the total tempera- 

ture 7. 


(mF me) d(mV + pA) 
a|——\ = - _ 
mv/ gRT mvV/¢gRT 
1 mV + pA dT 
2 mVgRT T 


(6 


The differential drag force dD, is given by the Fanning 
equation as 
dDy = F(1/2)pV°*A(4dx/D,) (7) 


where F is the friction factor; p, the mass density; x, 
the distance along the flow passage; and D,, the equiva- 
lent diameter of the passage. 

The differential heat-balance relation, which equates 
the heat transferred from wall to fluid to the heat ab- 
sorbed by fluid, is written in a form that permits direct 
substitution of the friction factor for the heat-transfer 
coefficient through the use of Reynolds analogy. 


" ~ [A/copgV : rare X 
gi = 2 (" F/9 ) (T,, — T)Fd (*) 


where / is the heat-transfer coefficient; c,, the specific 
heat at constant pressure; and 7,, the passage-wall 
temperature in degrees absolute. 

From Eqs. (1) and (5) through (8) the variation of the 
total-momentum parameter during the heat-exchanger 
flow process is expressible as 


m J ue pA 


a x eed 1 mV + pA 


2T mvWV/¢eRT 





(9 


stations in the flow passage. Thus, any station in a flow 
passage can be considered as a possible entrance of some 
other passage, so that if we know the flow variations in 
one flow passage we actually have answers for a large 
number of other passages which are segments of that 
original passage. With this point of view, it is only 
necessary to calculate the flow variations in a relatively 
few generalized flow passages and to plot the results for 
convenient determination of the flow variations in any 
segment of the generalized flow passages. This proced- 
ure was used to construct charts for handling the spe- 
cific case of air heated at constant passage-wall tempera- 
ture under turbulent-flow conditions. The method of 
constructing the charts and the manner of using the 
charts for obtaining solution of a typical heat-exchanger 
flow problem are briefly described herein. 
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The assumptions used in the construction of the 
charts are as follows: 

(1) The grouping (//c,pgV)/(F/2) is taken constant 
in accordance with Reynolds analogy. For air, this 
constant is calculated as 1.186. 

(2) y is taken constant as 1.400. 

(3) Fis assumed invariant with distance along the 
flow passage. 

On the basis of these assumptions, Eq. (8) can be 
formally integrated for constant 7, to give: 


i nae a (1 — 7) [1 — e—2:372 Fx/D)) (10) 
Ts Te Ve 

Eq. (10) relates the air-to-wall temperature ratios at 

two stations in a flow passage to the value of F(x/D,) 

between the two stations. In Eq. (10) subscript 1 de- 

notes the upstream station and subscript 2 denotes the 

downstream station. 

Integration of Eq. (9) must be performed numeri- 
cally; however, the integration is required for only a 
few generalized flow passages. The choice of the en- 
trance conditions for the generalized flow passages is 
such that all the possible flow and heat-transfer condi- 
tions of interest (7,,,/T from 5.0 to 1.1 and Mach Num- 
bers from 0.20 to choke) are encountered somewhere in 
the generalized passages. 

Although a single chart would have sufficed, the re- 
sults are presented in several working charts (Valerino') 
to obtain increased plotting and reading accuracy. The 
method of presentation of charts is illustrated in Fig. 1, 
which covers the following range of variables: (a) 
ratios of wall temperature to air temperature (degrees 
absolute) from 5.0 to 2.0, and (b) total-momentum 
parameters equivalent to Mach Numbers from 0.15 to 
1.00 (choke). The results are plotted relative to the 
conditions at a reference station herein defined as that 
station in the flow system at which the ratio 7;,/T is 
equal to 5.0. The reference station and all conditions 
at the reference station are denoted by the subscript 0. 
The ordinate on the chart is the ratio of the total 
momentum at any station in the flow system to the 
total momentum at the reference station. The abscissa 
is the reciprocal of the total-momentum parameter at 
the reference station. The curves are for different fric- 
tion length-diameter ratios F(x/D,), where x is the 
distance along the flow passage as measured from the 
reference station. The curve representing the choke 
limit for the flow is indicated on the chart. 

Toadd physical significance to the procedure involved 
in the construction and the use of the charts, a sketch 
of a generalized flow passage is given by the dashed lines 
at the bottom of Fig. 1. The position of the reference 
station, which as previously stated is characterized by 
T,,/T = 5.0, is indicated in the sketch. Different 
abscissa values on the chart correspond to different flow 
conditions at the reference station in the generalized 


passage. 
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Fic. 1. Pressure-drop chart for range 7y/7' = 5.0 to 2.0. 


Given the flow conditions at the reference station, the 
abscissa value on the chart is constant, so that as we 
progress from the reference station along the generalized 
passage we travel vertically downward from the curve 
for F(x/D,) = 0 on the chart. It is evident, then, that 
the locus of all flow conditions in any given flow passage 
lies on a common vertical line on the chart. Changing 
the flow conditions at the reference station merely shifts 
the abscissa value of the vertical line. The ordinate 
value, corresponding to the intersection of the vertical 
line and an F(x/D,) curve, gives the ratio of the total 
momentum at the station located F(x/D,) downstream 
of the reference station to the total momentum at the 
reference station. 

As F(x/D,) is increased, the total temperature 7 of the 
fluid is also increased, and thus the ratio 7,,/T is de- 
creased. Inasmuch as a unique relation exists between 
T,,/T and F(x/D,), each curve of constant F(x/D,) on 
the chart also represents the locus of points at which 
T,,/T is constant. Because the values of 7),/7 cor- 
responding to the values of F(x/D,) can be easily calcu- 
lated from Eq. (10), only the end values of 7;,/T are 
indicated on the chart; for F(x/D,) = 0.20, T,/T = 
2.0, and for F(x/D,) = 0, we are at the reference station 
so that 7,,/7 = 5.0. 

In a typical heat-exchanger problem, the flow condi- 
tions at the entrance of a passage of known length and 
diameter and specified wall temperature are given in 
which case it is required to determine the flow conditions 
at the passage exit (in particular, the static and total 
pressures and the total temperature). The operations 
involved in obtaining solution of a typical heat- 
exchanger problem by means of the charts are illustrated 
as follows: 

The given flow passage is considered as a segment of 
a generalized flow passage as indicated by the full lines 
in the sketch. The entrance and exit of the given flow 
passage are indicated in the sketch as stations en and 
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ex, respectively. The friction length diameter ratio of 
the given passage is F(L/D,), where F is the average 
friction factor for the passage as calculated for the given 
flow conditions, LZ is the passage length, and D, is the 
passage equivalent diameter. 

The first step is to locate the point on the charts cor- 
responding to the entrance flow conditions of the given 
passage. In essence, this step involves determination 
of, first, the value of F(x/D,) between the reference 
station and the passage entrance {referred to as [F- 
(x/D,)]en in Fig. 1} and, second, the value of the total- 
momentum parameter at the reference station and thus 
the chart abscissa such that the flow conditions obtained 
at the entrance of the given flow passage match those 
given in the problem. 

The value of [F(x/D,)]-, is calculated from Eq. (10) 
as that value of F(x/D,) required to effect a change in 
T,,/T from the reference value of 5.0 to the value speci- 
fied for the entrance of the given flow passage. The 
chart abscissa is then established by a simple construc- 
tion procedure as illustrated on the chart of Fig. 1. A 
straight line is drawn through the origin of coordinates 
with a slope that is easily calculated from the known 
entrance conditions of the given flow passage. The 
abscissa corresponding to the intersection of this slant 
line with the curve for F(x/D,) = [F(x/D,) |e, is the 
proper chart abscissa for the conditions of the problem. 
The point of intersection is the point on the chart repre- 
senting the entrance flow conditions of the given pas- 
sage. The ordinate for this entrance point is (mV + 
PA) en/(mV + pA)o, the ratio of the total momentum 
at the flow-passage entrance to the total momentum at 
the reference station. 

The next step is to locate the point on the charts 
corresponding to the exit flow conditions. This point 
obviously lies on a vertical line drawn through the 
entrance point—that is, it is located at the same chart 
abscissa as the entrance point. As is evident from the 
sketch at the bottom of Fig. 1, the F(x/D,) value for the 
exit point is simply equal to the sum of [F(x/D,)]en 
and the F(L/D,) of the given flow passage. Thus, the 
location of the exit point is fixed as indicated on the 
chart. 

The ordinate for the exit point is (mV + pA), + 
(mV + pA)o, the ratio of the total momentum at the 
flow-passage exit to the total momentum at the refer- 
ence station. The ratio of the ordinates for the exit 
and entrance points is thus the ratio of the total mo- 
mentum at the flow passage exit to the total momentum 
at the flow-passage entrance. 

Finally, we require an interpretation of the chart 
results in terms of such fluid-flow conditions as static 
and total pressures, velocity, Mach Number, etc. 
From the known entrance conditions, temperature rise 
as calculated by Eq. (10), and total-momentum ratio 
determined from the chart, the value of the total- 
momentum parameter at the flow-passage exit can be 
calculated. It is then a simple matter to obtain all the 


flow conditions at the flow-passage exit from suitable 
plots of the basis flow relations given by Eqs. (2) 
through (5). 

If F(x/D,) for the exit point is greater than 0.20 (the 
maximum value on the chart), the exit point appears on 
another chart, in which case it is necessary to transfer 
to the other chart in obtaining solution of a problem. 
Means for conveniently making this transfer have been 
incorporated into the charts. 

Heat-transfer and pressure-drop data have recently 
been obtained at N.A.C.A.* in an investigation con- 
ducted with air flowing through a smooth tube that 
was electrically heated to average wall temperatures of 
from 50° to 1,200°R. above the entrance-air tempera- 
ture. These data show an important effect of tempera- 
ture differential between tube wall and fluid on the 
heat-transfer and friction phenomena. The effect on 
heat transfer has been accurately isolated for Reynolds 
Numbers above 10,000 as expressed by the correlation 


equation 


AD, Ved. 08 . . 0.4 
> = 0.023 (* £ ) (‘xt (11) 
Ry Kw kw 


where subscript w denotes fluid conditions evaluated at 
the average tube-wall temperature, k is the thermal 
conductivity, and yu is the absolute viscosity. The other 
symbols are as previously defined. The significance of 
pwV in the modified Reynolds Number term p, VgD,/ 1. 
is given by 

PwV _ pV (pw/p) = eV(T Ty) (12) 


In other words, p,. |’ is the conventional bulk velocity 
pV multiplied by the ratio of the fluid density at the tube 
wall to the fluid density of the bulk stream p,,/p; the 
assumption is then made that p is inversely proportional 
to T rather than ¢. 

Comparison of Eq. (11) with the standard heat- 
transfer correlation equation presented by McAdams’ 


h 7 (F=\" (a=) (t= )" 13) 

hs Be Cool Win cg 
where h/h, is the ratio of the actual to the standard 
heat-transfer coefficient. The Prandtl Number for 
air was assumed constant in obtaining Eq. (13). The 
subscript av denotes fluid conditions evaluated at the 
average fluid temperature within the flow passage; the 
average values of 7, c,, and u used in Eq. (13) are con- 
sistent with the determination of an average value of / 
in the flow passage. 

From a partial analysis of the pressure-drop data 
obtained in the experiments, Reynolds analogy appears 
valid in the range of high-temperature differential be- 
tween wall and fluid. It is accordingly assumed that 
Reynolds analogy applies—that is, that a change in the 
heat-transfer coefficient is accompanied by a propor- 
tionate change in the friction factor. This assumption 
is expressed in Eq. (14): 
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Entrance Conditions 


Average 
Tw, : 
R. Reynolds P, |bs. per a» Mach 
Number sq.ft. “ih. Number 

1,672 | 184,000 4,527 528 0.37 
1,072 155,000 3,682 554 0.41 

861 219,000 4,458 540 0.46 
1,670 263,000 6,186 530 0.39 
1,752 | 53,600 2.417 527 0.21 


Test Calculated Percentage 
Exit Results Results | Difference 
Mach ————— a in 

Number Pex Ad, lbs. Pez Ap, lbs. 
pen | per sq.ft. Pen per sq.ft. Ap 

0.84 0.551 2,031 0.544 2,065 ® 3 
0.70 0.639 1,329 0.657 1,264 —4.9 
0.87 0.568 1,926 0.553 1,993 3a |] 
1.00 0.464 3,315 0.458 3,353 1.29 
choke 


0.32 0.881 288 0.889 267 —7.1 


Calculated choking length less than actual choking length by 2.6 per cent. 


Fic. 2. Comparison of calculated and test results. 


F/F, = h/h, (14) 


where F/F, is the ratio of the actual friction factor to 
the standard friction factor used in current heat- 
exchanger practice. Based on this assumption, a com- 
parison was made, for several experimental runs, be- 
tween the pressure drops measured in the tests and the 
pressure drops obtained by use of the charts for the 
entrance flow and average wall-temperature conditions 
of the experiments. This comparison is tabulated in 
Fig. 2. The range of conditions covered in the compari- 
son are: (a) tube-wall temperatures from 861° to 
1,752°R.;  (b) Reynolds Numbers from 53,600 to 
263,000; (c) Mach Numbers at tube entrance from 
0.21 to 0.46; and (d) Mach Numbers at tube exit from 
0.32 to 1.00 (choke). The static pressure ratios across 
the tube are from 0.881 to 0.464, the static pressure 
drops from 288 to 3,315 lbs. per sq.ft. The percentage 
differences between the calculated and the test pressure- 
drop values vary from 3.5 to —7.1 per cent; this agree- 
ment is considered good. At the relatively low Rey- 
nolds Number of 53,600 the agreement is not so close 
as at the higher Reynolds Numbers. This result is 
attributed to the laminar-flow boundary-layer régime 
at the tube entrance, which, because of the smooth well- 
rounded entry of the test tube, occupies an extensive 
portion of the tube length even for relatively high values 
of Reynolds Number. The good agreement obtained 
between the experimental and the chart pressure-drop 
values indicates that the assumptions of Reynolds 


analogy and of one-dimensional flow (used in the 
basic relations) are valid even for the high subsonic flow 
Mach Numbers and high rates of heat input to the 
fluid. 

In summary, an analysis has been made of the com- 
pressible flow variations occurring in heat-exchanger 
passages. Based on this analysis, working charts for 
the direct determination of pressure drop have been 
constructed for the specific case of air flowing turbu- 
lently through smooth constant-area passages wherein 
heat is added to the air stream by the passage walls that 
are assumed to be at a constant temperature throughout 
their length. The accuracy of the charts and assump- 
tions involved therein have been verified by checks 
with experimental data for extreme conditions of high 
subsonic Mach Number and high rates of heat input 
to the fluid. 
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A Hypothesis Concerning Turbulent Exchange 


Wallace D. Hayes 
Assistant Professor, Graduate Division of Applied Mathematics, 
Brown University, Providence, R.1. 


January 31, 1949 


— EXCHANGE PROPERTIES of an incompressible fluid in a 
state of anisotropic turbulence are characterized by the 
symmetric exchange tensor A;;._ This quantity is the generaliza- 
tion of the exchange coefficient of isotropic turbulence. With re- 
spect to principal axes this tensor is expressed by three compon- 
ents only, Ay, Ag, and A33, the off-diagonal terms being null. 
These components are all positive. The hypothesis is that the 
inequality 
|Au ~ Az | < Ass 

holds, where, of course, the indices may be permuted. This rela- 
tion is in strict analogy with the similar relation for the moment 
of inertia tensor for rigid bodies. 

This hypothesis may be made to appear reasonable with a few 
qualitative considerations. The exchange tensor component Ay 
is generally expressed 


Ay - pul, 


where x; is the root-mean-square velocity fluctuation in the direc- 
tion , and J, is a suitably defined mixing length for the same 
direction. This view takes no account of the effect of the con- 
tinuity condition in connecting the motions in various directions. 
A better expression would be 
Ay = pa2R,? + pw3R;3? 

where w: is the root-mean-square vorticity of an eddy whose axis 
is in the direction » and R, is a suitably defined mixing radius 
describing the size of the eddy. The hypothesis given above is an 
immediate result of this expression, 
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Complementary Energy Procedure for Flutter 
Calculations 


Eric Reissner 

Associate Professor, Massachusetts Institute of Technology, 
Cambridge, Mass. 

February 1, 1949 


IX THE FOLLOWING there is outlined a procedure for flutter calcu- 
lations which bears the same relation to the well-known 
Rayleigh-Ritz procedure! as does the theorem of minimum com- 

1 See, for instance, U.S. Air Forces Technical Report No. 4798 (1942) by 
B. Smilg and L. S. Wasserman. 


plementary energy in elasticity (minimum principle for the 
stresses) to the theorem of minimum potential energy (minimum 
principle for the displacements). 

The procedure is stated here for the case of torsion-bending 
flutter of a cantilever wing. The possibility of its extension to 
more general flutter problems will be apparent. 

The differential equations for torsion-bending flutter of a canti- 
jever without sweep may be written in the form 


d?*M dx? = —w* (any + 426) (1) 
dT /dx = —w* (aay - A220) 2) 
M = — El d*y/dx?, TT = GJ d0/dx (3 


where the coefficients a,,,., depend in a known way! on reduced 
frequency and other parameters. 
Boundary conditions may be taken in the form 


x = 0, y = dy/dx = 0=0 4 

x=L, M =dM/dx = T =0 5) 

The variational principle that forms the basis of the proposed 

procedure may be stated in the following way: Among all states 

of stress (M, T) and displacements (y, 0) that satisfy Eqs. (1), (2), 

and (5), the state that also satisfies Eqs. (3) and (4) is determined 
by the variational equation, 


p[ MoM , Ter 
Le * @ 


— w*y(audy + 4260) — 


w?O(andy + a) [as = 0 (6) 


The proof of this theorem consists in showing that Eq. (6) is 
satisfied if one substitutes in it Eqs. (3) and takes account of the 
fact that 6./, 67, dy, and 56 must satisfy Eqs. (1), (2), and (5), 
as well as V, T, y, and @. 

To use the theorem for the approximate determination of flut- 
ter speeds, set 


y = AY, 6 = Bo (7) 


where Y and 0 are assumed modes—for instance, the fundamental 
uncoupled bending and torsion modes of the cantilever beam or, 
more simply, Y = (x/L)? and 6 = (x/L)—and where A and B 
are constant amplitude factors. Substitute Eqs. (7) in Eqs. (1) 
and (2) and calculate 
M = w3(AM, + BMo)( g) 
T = w?(AT, + BT) ) vee 
such that M and T satisfy Eqs. (5) for all values of A and B. 
Corresponding to Eqs. (7) and (8) one has dy = (6A) Y,etc. The 
results are substituted in Eq. (6), and this equation reduces them 
to the following two simultaneous equations for A and B 
kyA +. kw»B = Ol 
ky A aa kB = 0) 


The coefficients k are given by the following expressions 


Jo : E (= + 72) — an r] dx (10a 
ky» = Ey [. (Seve = “= ) — ay ro | dx (10b 
kn = i E (4% “t- at) — on | dx (10c 
kn = fy” [ (Ze + 7) ~ on0* lex (10d) 


The determinantal equation determining the flutter speed 
follows from Eqs. (9) to 


kiko — Rvkn = O (11) 


Ru 


It is believed that the use of Eqs. (10) and (11) may lead, for 
the same modes Y and 9, to more accurate results than the Ray- 
leigh-Ritz procedure, because in the present procedure the strain 
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energy is expressed in terms of functions that are obtained from 
the assumed modes by integration rather than by differentiation 
as in the Rayleigh-Ritz procedure. This may permit the working 
with less close approximations to the actual modes than is pos- 
sible with the Rayleigh-Ritz procedure; or it may, in some in- 
stances, permit to take account of fewer degrees of freedom than 
is necessary with the Rayleigh-Ritz procedure. These advantages 
will be somewhat offset by the fact that evaluation of the coef- 
ficients k involves a larger number of integrations than does 
evaluation of the corresponding coefficients in the Rayleigh-Ritz 


procedure, 


Analysis of an Unsymmetrical Pin Pattern 


R. D. Hiscocks 
The de Havilland Aircraft of Canada, Ltd., Toronto, Canada 


February 16, 1949 


I THE PAPER “Analysis of Eccentrically Loaded Riveted or 

Bolted Joints,’’ presented by G. M. Monroe in the February, 
1949, JOURNAL OF THE AERONAUTICAL SCIENCES, attention is 
drawn to a useful method of locating the most highly loaded pin 
in a rivet or bolt cluster. A variant of this method will reduce 
even further the amount of arithmetic required to analyze an 
unsymmetrical pin pattern. 

With reference to Fig. 1, the procedure is as follows: 

(1) Locate G, the center of resistance of the cluster on the 
basis of effective bearing or shear areas. This can most readily be 
accomplished graphically on a full-scale drawing. (See Basic 
Structures by F. R. Shanley.) 

(2) Extend the line of action of the applied load P, and from 
this line erect a perpendicular that passes through the centroid 
G and extends on the side of G remote from P a distance GQ, so 
that 

GQ = SAr?/eTA 
where 

A effective area of pin in shear or bearing, 

r = radial distance of pin from center of resistance, and 
the distance from line of action P to center of resistance. 


e 
(8) Determine the radial distance L of any pin from Q. The 
load in that pin is then 
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The most highly stressed pin is obviously that for which L is a 


PeAL/ZAr* 


maximum. 

Proof of the method, which I believe originated with Prof 
W. G. Sutton (Journal of the Institution of Structural Engineers, 
September, 1933) follows from elementary geometry. 

In calculating the centroid, the question is often raised whether 
to use effective bearing or shear areas. With the ratio of rivet 
diameter to shear thickness employed in Mr. Monroe’s example, 
it would appear preferable to use bearing areas, since the indi- 
vidual rivet displacements at limit load would be predominantly 
in bearing. 


Heat Addition to a Flowing Gas 


Newman A. Hall 

Professor of Thermodynamics, Mechanical Engineering Department, 
University of Minnesota, Minneapolis, Minn. 

February 16, 1949 


Is THE DISCUSSION “‘On the Addition of Heat to a Gas Flowing 

in a Pipe at Subsonic Speed,”’ by Joseph V. Foa and George 
Rudinger (Journal of the Aeronautical Sciences, February, 1949) 
the authors have strongly implied that extended revision is in 
order for most of the previous material appearing on this subject. 
The extensive successful experimental confirmation of these previ- 
ous theories suggests further examination. 

As pointed out in the paper, it is important in compressible 
flow to keep in mind the history of the establishment of the par- 
ticular flow configuration under consideration. However, once 
established, steady flow will obey laws that have been stated 
correctly elsewhere.! 

It appears that the authors have leaned too far in the other 
direction in some of their inferences. For example, they state 
that, when flow is initially sonic, there may be an upstream static 
temperature rise in excess of the downstream rise because of the 
addition of heat. The inference follows that the upstream tem- 
perature rise is due to a transmission of heat upstream. This, of 
course, is not true. Under the usual assumption of neglecting 
thermal conductivity, which the authors themselves tacitly 
make, heat added cannot move upstream of the streamline, in a 
time-distance plane, where it is added. Further, an identical 
upstream static temperature rise could be obtained by introduc- 
ing a suitable downstream obstruction. It would be more proper 
to say that the upstream static temperature rise was due to the 
flow resistance associated with heat addition. 

The implication is made that, in talking of heat addition, the 
phrase ‘‘initial conditions” must refer only to the state in time 
preceding addition of heat anywhere in the system. Since many 
investigators are concerned primarily with steady state condi- 
tions where heat is being added at a constant rate, it has been 
convenient to refer to the state in space preceding heat addition 
as the initial condition. This is a valid use of the term and, 
insofar as a given element of fluid is concerned, represents initial 
conditions in time as well. Accepting this definition, which is 
more widely established, rather than that of the authors, certain 
conclusions cannot be questioned. In particular, the following 
statements remain true: 

(1) In steady one-dimensional flow, when heat is added at a 
constant rate, the Mach Number of the flow will increase steadily 
from its initial value through the region wherein heat is added 
but, in no case, can exceed a value of 1 if the flow is initially sub- 


sonic. 


' Shapiro, A. H., & Hawthorne, W. R., The Mechanics, and Thermo- 
dynamics of Steady One-Dimensional Gas Flow, Journ., Appl. Mech., Vol. 
14, pp. 317-336, December, 1947. 
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2) In steady one-dimensional flow, when heat is added at a 
constant rate, the static temperature will increase through the 
region where heat is added, provided the initial Mach Number, 
M, < 1/</y, but will decrease when the initial Mach Number 
lies in the interval, 1/./y < M <1. 

Both of these conclusions have been correctly stated and devel- 
oped in many places—e.g., reference 1—and are also clearly and 
unavoidably contained in the analysis presented by the authors 
themselves. 

It has been long accepted terminology to describe a flow in a 
convergent divergent nozzle as choked when the throat Mach 
Number is 1. The only restriction on flow thus implied is that a 
change in downstream conditions cannot increase the mass flow. 
Furthermore, in the commonly accepted sense the choked condi- 
tion represents the maximum mass flow consistent with fixed 
initial conditions. Conversely, the size of a choked nozzle cannot 
be reduced without modifying the initial conditions. The choking 
action of heat addition parallels in every respect that of flow re- 
striction in a nozzle. 


i 


Errata 


Herman F. Michielsen 
Fokker Aircraft Factories, Amsterdam, Holland 
February 7, 1949 


— FOLLOWING CHANGES are submitted for my article ‘‘The 
Behavior of Thin Cylindrical Shells After Buckling Under 
Axial Compression”’ (Journal of the Aeronautical Sciences, Vol. 
15, No. 12, pp. 738-744, December, 1948). 

Page 739—The second line of Eq. (2) should read 


2mx 4 2ny (2) 
7g | COS ; 2 
go | C R cos R 


Page 739—The second and third lines of Eq. (3) should read: 


y1+ 16u4 dys 
n' 14 + 
's G+? at oe? 


4 ri | 
= ; ] £17g27 + = mm. aE +... 
(9 + 2)? (1 + u2)2?™ f 2 
The same error in Eq. (3) appears in Eqs. (5), (6), (8), (9), 
(12), (18), (28), and (29); in the two equations between Eqs. (28) 
and (29); and, finally, in Eqs. (30) and the one immediately above 
Eqs. (30). 


Page 741—The equation at the top of the page should read: 


oe aes ti" tm)... 
27 4 


Eq. (17) should read: 


. ih = Be 
(1 + p*)*y (« +yv+ { = : ta.) (17) 


Eq. (21) refers to the latter of the two equations connected 
with a brace. 


Page 742—The first line of Eq. (29) should read: 


1 4y Sw? 
1 . - 
we opt(l + y*) ¥ 


Page 743—The last equation on the page should read: 
Tz,y,b = — Et/(1 + v)(0*w/dx dy) 
Page 744—The first of Eqs. (31) should read 


oz, oR 


Ei = -..[(1 +2) co... 
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The Shock-Boundary Condition 


B. Etkin 
Assistant Professor of Aeronautics, University of Toronto, Canada 
February 28, 1949 


POINT THAT APPEARS to have been overlooked in the litera- 
ture on supersonic flow is that a condition exists which 
makes possible the solution for the flow aft of a curved shock and 
at the same time serves to determine the form of the shock itself, 
Consider the situation illustrated in Fig. 1. Two states are 
separated locally by a shock wave, the solution in state 2 being a 
function of that in state 1 and the shock parameter a. This 
statement may be expressed by 


(state 2) = f [ (state 1), a] 1) 





SHOCK 







STATE 2 


Fic. 1. 


In general, it is possible to construct a characteristic solution 
for state 2, which grows outwards from the boundary. Let it be 
assumed that the solution is known at S. Hence the charac- 
teristic segment C,° is known, and also the characteristic condi- 
tion along it. This characteristic condition is a relation among 
the variables of state 2, which we write 


f (state 2) = 0 2 


When state 1 is given, then Eqs. (1) and (2) determine a and 
state 2. If the technique indicated here be applied stepwise 
from the point of origin of the shock wave outwards into the flow 
region, the whole solution may be constructed. 

The characteristic equations for state 2 may be used in the 
nonisentropic or isentropic form, depending on the precision re 
quired. The solution at the point P will normally have to be 
carried out numerically or graphically, since the equation that 
results from substitution of Eq. (1) into Eq. (2) is extremely 
complicated. 

When applied to the case of a circular are airfoil of thickness 
ratio 10'/s per cent, the solution for the pressure differed from the 
linear theory by 17 per cent at the loading edge and gave a drag 
5 per cent greater. 





Note on the Perturbation Theory for Axially 
Symmetric Flows 


C. C. Lin and E. Reissner 
Massachusetts Institute of Technology, Cambridge, Mass. 
February 3, 1949 


iw PRESENT NOTE is concerned with the unsteady motion in 
the transonic range of slender bodies of revolution. An 
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analysis has been undertaken along the lines of some earlier work 
on the problem of two-dimensional flow.! 

The following question appears to be of particular interest: 
Under what conditions is it possible to treat the problem of ac- 
celeration through the transonic range by means of a linearized 
theory? We have found that the magnitude of the acceleration a 
necessary to permit application of a linearized theory must satisfy 
the following order of magnitude relation, 


a/g > (co?/gb)dp (1) 


In Eq. (1) the quantity co is the velocity of sound in the undis- 
turbed medium, g is the acceleration of gravity (taken arbitrarily 
as reference), b is a typical linear dimension of the body along its 
axis, p is the ratio of a typical cross sectional radius of the body 
to b, and 6 is the ratio of a typical waviness of the body to b. 
For a slender body in the usual sense one has p = 6 and then the 
condition reads a/g >> (cy?/gb)&? 


The range of values of the Mach Number .V/ in which the above 
conditions are required for linearization is given by 


1— M)| > O(ép) 
and by 

1— M| > O(8) 
respectively. As soon as the values of .V are sufficiently different 
from unity so as to lie outside these ranges the problem can be 
linearized regardless of the magnitude of the acceleration. 

Eq. (2) indicates for instance that, when the length 6} of the 
body is 10 ft. and the slenderness ratio 5 has the value 0.1, then 
the necessary acceleration a must be large compared with 30g. 
It may be seen from the order of magnitude relation involving 
M that the range of / which is of importance here is given by 

1 — M| = O(0.01) and is thus quite narrow. 

Derivations of the above results and of other results for axially 

symmetric and two-dimensional flows will be reported on in the 


near future. 


1 Lin, C. C., Reissner, E., and Tsien, H. S., On Two-Dimensional Non™ 
Steady Motion of a Slender Body in a Compressible Fluid, Journal of Mathe~ 
matics and Physics, Vol. 27, pp. 220-271, 1948 


End-Fixity Coefficients for Columns 


Joseph S. Newell 
Department of Aeronautical Engineering, M. |. T., Cambridge, 
Mass. 


February 25, 1949 


HE FEBRUARY ISSUE OF THE JOURNAL OF THE AERONAUTICAL 

SCIENCES contains two papers pertaining to the estimation 
of end-fixity coefficients for columns. A relation developed from 
Equation 14:6, page 67, Vol. II of Airplane Structures leads to 
similar results with, perhaps a little less computational labor. If 
it be assumed that a practical column will deflect before it is sub- 
jected to its critical end load and that such deflection entails the 
development of moments in each end supporting structure which 
tend to oppose increasing deflections, the beam-column side load 
w in equation 14:6 may be set equal to zero to give an expression 
for the deflection y at any section x in. from the origin. It is 


I Mz, — M, M, — M, cos (L/j) 
y=5| + =—)s- cos (4/5) 


sin (L/j) 
If the shape of the deflected column be such that the distance 
between points of inflection is L/n in., the critical load on the 
column will be P,., = n?x*(EI/L?) so long as all stresses are in the 


sin (x/j) — M, cos (x ‘| 
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elastic range. Substitution of P,, for P in the quantity 7 = 
V EI/P yields nw for L/j and nx(x/L) for x/j. 

Since dy/dx = a, the slope of the deflected column, whereas 
d*y/dx? = M/EI, it is a simple operation to obtain the ratio 
M/EIea and eliminate the axial load P or P,, from the expression. 
If this ratio be evaluated for M = Mh, a = m, at x = Oit will be 
found that 4, factors from the numerators on both sides of the 
resulting equation so the expression yields 1/EJay. 

Substitution of MW = Ms, a = az, at x = L yields a similar ex- 
pression for 1/EJa2. Elementary manipulation of the reciprocals 
of these expressions gives: 

nr sin nr z (* + ) 


= = ¢(N) 
(cos nr — 1) El , 


a2e2—~— ai 


Computations of g(N) for various values of n gives: 


n 0.1 0.3 0.5 0.7 0.9 

¢(N) —1.98 —-1.85 -—-1.57 -— 1.12 —0.448 
n ‘3 1.3 1.5 Re 1.9 
¢(N) +0.547 +2.08 +4.71 +10.48 +37.70 


The quantities L, EZ, and J have their usual significance and per- 
tain to the member under investigation. If that member be taken 
to have its longitudinal axis horizontal and if the signs of quanti- 
ties comply with the beam-column conventions, deflection y is 
positive upward; M, is positive when clockwise at the left end; 
Mz. is positive when counterclockwise at the right end; a; and ae 
are positive if the axis of the column rotates counterclockwise 
from its unloaded position, so a positive .4, produces a negative 
a and a positive ./2 produces a positive as. 

The end moments J/, and J, may have any convenient values 
when computing a; and az, since it is the ratios that are involved in 
¢g(N). It is important to note that the column is assumed to be 
carrying its critical load so it is neutrally stable and has no ability 
to transmit external forces or moments. Its ends assume the 
slopes taken by the end supporting structures. The shape of the 
deflection curve depends in part on those slopes and in part on 
the effects of \4,, Mo, and P, at each section. So long as the mem- 
bers composing the supporting structure carry no axial loads and 
have their ‘‘far-ends’”’ fixed, their slope changes may be found from 
a = M,>OL/4EI and a = M:>°L/4EI. If the members do 
carry axial loads, the stiffness factors L/4E/ must be corrected for 
the increases or decreases resulting from such loads, and if the 
far-ends are not fixed, the factor 4 must be modified. 

For the bent shown in Fig. 3 of G. C. Best’s paper (JOURNAL OF 
THE AERONAUTICAL SCIENCES, Vol. 16, No. 2, p. 117) the axial 
loads in members 1-3 and 2—4 are zero and the far-end of each is 
fixed, so the stiffnesses are L/4EI. 

For M, = 104, 


M,Li-; 10 X 50 
a= =— = —3 13 
4(E7), 16 X 104 
and 
10* X 50 4.1.95 
a = + 4x10 <0 
For column 1-2 the factor L/EI = 30/10‘ so: 
30 (104 + 104 60 X 10% 
¢(N) = - — = - = 13.7 
10% 1.25 + 3.13 4.38 X 104 
From a plot of ¢(.V) versus n, nm = 1.755, which corresponds with 


Mr. Best’s Ve as it should, since the restraint coefficient c is 
equal to the square of the number of half-waves in the deflection 
curve. 

Negative values for ¢(N) indicate that 44, and M, tend to 
increase the deflection of the member rather than to oppose it 
and that the shape of the deflection curve is such that the effective 
pin-ended length L./n is greater than the actuallength LZ. Such 
values are possible, though they appear infrequently. The 
standard case that comes to mind is the ‘‘axially-loaded canti- 
lever,’’ the column with one end free, one fixed, for which n = 
1/2,c = 1/4, e¢(N) = —1.57. 
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